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aO  abstract  f C tut  tint  ic  on  reverie  eide  II  necenary  mrd  Identity  by  block  number) 

- A finite  element  program  is  described  for  the  analysis  of  subsonic  and 
transonic  flows  over  thin  airfoils  by  solving  the  nonlinear  transonic  poten- 
tial equation  based  on  small  $1 j,sjturbance  theory.  The  present  numerical  al- 
gorithm Va  developed uatri^me  concept  of  finite  elements  in  conjunction  with 
the  least  square  method  of  weighted  residuals.  Since  the  governing  equation 
is  of  the  elliptic -hyperbolic  type,  a I’one-sided  assembly  technique*  was  de- 
vised and  adapted  in  the  supersonic  region  to  restore  the  directional  pro- 
perty of  the  flow,  which  was  removed  by  the  exclusion  of  entropy  from  the 
transonic  potential  equation.  The  finite  element  discretization  results  in  a 
system  of  banded  nonlinear  algebraic  equations,  which  is  solved  by  direct 
iterations.  The  elements  presently  used  include  triangles  and  quadrilaterals 
with  the  perturbed  potential  function  and  velocity  components  as  nodal  un- 
knowns. Boundary  conditions  of  both  Dirichlet  and  Neumann  types  are  there- 
fore imposed  convex.'  ently.  'Alsn/^secondary  unknownst  such  as  local  Mach 
number  and  pressure  coefficient,  etc.  f are  computed  directly  without  resort- 
ing to  numerical  differentiation.  The  ^pr-esanlfcomputer  program  is  separated 
into  two  parts:  the  first  part  (designated  as  STR AN L-J)  generates^  from  a 
limited  number  of  input  cards, cTKe  necessary  mesh  information  and,  if  desired, 
produces  a mesh  plot  and  optimal  nodal  numbering  as  well;  the  second  part 
(STRANL-II^Carries  out  the  analysis  and  displays  the  pressure  coefficients 
along  the  chord  line  on  printer  plots.  Two  sample  cases  of  flow  over  a NACA- 
64  A006  and  a 6%  thick  circular  arc  are  given^to  demonstrate  the  applicability 
and  usage  of  the  program.  The  solution  procedures  are  found  to  be  very  ac- 
curate and  quite  efficient,  permitting  the  aerodynamic  forces  to  be  calculated 
to  engineering  accuracy  in  tens  of  seconds  on  a CDC  6600  computer. 
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1.  INTRODUCTION 

The  aim  of  the  present  study  was  to  develop  an  efficient  and  accurate 
numerical  algorithm  for  the  analysis  of  steady  transonic  flow  over  nonlifting 
thin  airfoils.  With  this  objective,  the  present  formulation  is  therefore  based 
on  the  small  disturbance  but  nonlinear  transonic  potential  equation  for  in  - 
viscid,  compressible  fluid. 


It  is  well  known  that  two  major  difficulties  exist  in  solving  the  transonic 
flow  problem:  (1)  the  governing  partial  differential  equation  is  nonlinear, 

and  (2)  the  equation  is  of  mixed  elliptic -hyp erbolic  type.  Because  of  these 
difficulties,  numerical  procedures  must  therefore  be  used,  with  provisions 
made  to  properly  account  for  the  above  transonic  effects.  Among  the  numer- 
ical approaches,  the  finite  difference  relaxation  technique  is  the  most  well 
developed  and  widely  used  in  recent  transonic  calculations.  On  the  other 
hand,  the  finite  element  technique,  in  spite  of  its  tremendous  success  in  the 
field  of  solid  and  structural  mechanics,  has  never  been  attempted  until  the 
present  investigation,  obviously  because  conventional  finite  element  technique 
is  invalid  in  solving  mixed-type  equations.  This  difficulty,  nevertheless,  was 
successfully  resolved  by  a special  assembling  technique  devised  for  the  super- 
sonic region,  while  the  major  advantages  of  the  conventional  finite  element 
method  are  retained.  These  include  the  complete  flexibility  in  element  ar- 
rangement, its  effectiveness  to  cope  with  complex  geometric  shapes  and 
boundary  conditions,  and  the  ease  with  using  efficient  higher  order  approxi- 
mations. For  these  reasons,  the  present  approach  appears  to  be  superior 
in  solution  accuracy  and  computational  efficiency  compared  to  other  exm-ting 
numerical  techniques. 

The  present  numerical  algorithm  is  developed  using  the  concept  of  finite 
elements  in  conjunction  with  the  least  squares  method  of  weighted  residuals 
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The  basic  element  presently  used  has  a cubic  expansion  for  the  per- 
turbed velocity  potential  inside  the  element,  with  nodal  unknowns  repre- 
senting the  function  itself  and  the  two  perturbed  velocity  components.  The 
resulting  system  of  nonlinear  algebraic  equations  is  banded,  although  non- 
symmetric.  This  system  oi  algebraic  equations  is  solved  by  direct  iteiative 
procedures. 

The  theory  and  numerical  procedures,  together  with  the  usage  of  the 
computer  programs,  are  discussed  in  the  following  sections.  The  theory  in- 
cluding the  small  disturbance  potential  equation  with  its  corresponding 
boundary  conditions  and  the  related  secondary  unknowns  are  summarized 
in  Section  2.  The  numerical  procedures  employed  in  the  present  approach 
are  discussed  in  Section  3.  These  include  the  finite  element  formulation 
based  on  least  squares,  elements  used  in  the  computation,  special  considera- 
tions for  the  supersonic  region,  treatment  of  boundary  conditions  and  a de- 
scription of  the  iterative  procedures.  In  Section  4,  the  two  parts  of  the 
computer  code,  namely,  STRANL-I  and  STRANL-II  are  described  separately, 
in  the  aspects  of  scope  and  flow  chart,  description  of  variables,  subroutines 

used,  and  finally  input  and  output.  Two  sample  cases  are  given  in  Section  5 
to  demonstrate  how  to  use  these  computer  programs,  which  are  listed  in  the 
Appendixes. 


1 

i 


2.  THEORY  - ASSUMPTIONS  AND  BASIC  EQUATIONS 

The  objective  of  the  present  study  was  to  develop  an  efficient  and  ac- 
curate numerical  algorithm  for  the  analysis  of  steady  transonic  flow  over 
nonlifting  thin  airfoils.  With  this  objective  in  mind,  the  formulation  was 
therefore  based  on  the  small  disturbance  but  nonlinear  transonic  potential 
equation  for  inviscid,  compressible  fluid.  The  embedded  shock  is  assumed 
to  be  weak  and  boundary  layer  effects  are  neglected.  These  assumptions 
render  the  small  perturbation  theory  possible  and  lead  to  the  following 
mathematical  problem  to  be  solved. 


Differential  Equation 


1 - M 


(1  + » <i> 


+ <j> 

xx  , yy 


Boundary  Conditions 


( 1 f u>  - v 


on  the  airfoil 


V = 0 


at  infinity 


v = 0 on  line  of  symmetry  (4) 

where  <£  = perturbed  velocity  potential  function,  u,v  =perturbed  velocity  com- 
ponents in  the  x-  and  y-directions,  respectively,  M = freestream  Mach  number, 
y = ratio  of  specific  heats,  taken  to  be  1.4  for  air,  and  g = function  of  x defin- 
ing the  geometry  of  the  airfoil.  Equation  (1)  is  in  dimensionless  form  and  the 


. . 


x-axis  is  aligned  with  the  undisturbed  flow  direction.  The  dimensional  (with 
primes)  and  nondimensional  quantities  are  related  by 


and  <f>  = — 

c Y U c 


where  c and  are  the  characteristic  length  and  speed,  which  are  currently 
taken  as  the  chord  length  and  the  freestream  speed,  respectively. 


Figure  1 shows  the  flow  field,  together  with  the  steady  transonic  govern- 
ing equation  and  corresponding  boundary  conditions.  Because  only  nonlifting 
thin  airfoils  are  being  considered,  the  airfoils  are  necessarily  symmetric  and 
hence  only  half  of  the  flow  field  is  depicted  and  analyzed.  As  the  figure  shows, 
the  flow  tangency  condition  on  the  airfoil  is  retained  in  its  nonlinear  form,  as 
defined  by  Eq.  (2),  and  is  to  be  applied  along  the  actual  airfoil  surface. 

Another  approach  is  to  use  the  linearized  boundary  condition  and  apply  it 
along  the  airfoil  chordline,  as  commonly  adapted  in  most  existing  finite  differ- 
ent e codes.  Apparently,  use  of  the  nonlinear  form  is  more  accurate  when 
the  magnitude  of  u becomes  non -negligible.  On  the  other  hand,  use  of  the 
linearized  boundary  condition  could  enhance  solution  convergence  somewhat 
as  the  nonlinearity  now  exists  in  the  governing  equation  only.  Extensive 
studies  and  comparisons  are  certainly  required  in  order  to  assess  the  two 
approaches.  In  the  present  computer  code,  the  nonlinear  form  is  used  unless 
otherwise  specified  through  an  input  control  key,  in  which  case  the  linear  form 
will  be  used.  On  the  line  of  symmetry,  as  is  well  known,  the  Neumann  condi- 
tion of  no  norma]  clow  is  to  be  imposed.  For  the  far  field,  because  only  a finite 
domain  can  be  treated  in  the  analysis,  the  condition  of  no  disturbance  at  infinity 
is  assumed  to  be  valid  on  the  outer  boundary  which  is  usually  taken  as  a few 
chord  lengths  away  from  the  airfoil.  Numerical  experimentations  indicate  that 
a flow  region  with  H ^ 2c  and  L>4c  is  generally  required  to  yield  solution  with 
adequate  accuracy.  For  flow  with  higher  freestream  Mach  number,  the  flow 
region  should  be  extended  somewhat  to  account  for  effects  from  the  enlarging 
supersonic  pocket. 


For  the  problem  of  lifting  airfoils,  the  entire  flow  domain  mus!  be  con- 
sidered, In  this  case,  an  asymptotic  solution  including  at  least  the  circulation 
as  unknown  should  be  applied  on  the  finite  outer  boundary.  This  is  part  of  a 
current  study  conducted  for  AFFDL  and  will  be  implemented  later. 


Once  the  flowfield  solution  in  terms  of  the  perturbed  velocity  potential 
has  been  obtained,  all  secondary  unknowns  can  be  calculated  subsequently. 
These  include 
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In  the  above,  = 1,  the  normalized  freestiream  speed,  a = local  speed  of 


sound,  M = local  Mach  number,  V = the  total  velocity,  p = local  static  pressure, 


Pq  = stagnation  pressure,  and  = pressure  coefficient. 
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3.  NUMERICAL  SOLUTION  PROCEDURES 

The  concept  of  finite  elements  in  conjunction  with  the  method  of  weighted 
residuals  (MWR)  is  the  basis  of  the  present  numerical  procedures  in  solving 
the  small  perturbation  transonic  flow  equation.  Of  all  the  MWR  methods,  the  I 

Galerkin  and  the  least  squares  approaches  were  investigated  (Chan  et  al., 

Ref.  1),  of  which  only  the  latter  approach  was  found  applicable  to  mixed  flow 
analyses.  The  success  with  the  least  squares  formulation  is  attributable  to 
the  resulting  matrix  being  positive  definite  and  well  conditioned.  This  approach 
was  therefore  incorporated  in  the  present  computer  program,  with  numerical 
procedures  described  in  the  following  subsections.  These  include  the  finite 
element  formulation,  element  description,  special  considerations  in  the  super- 
sonic region,  the  imposition  of  boundary  conditions  and  iterative  procedures. 

3.1  Finite  Element  Formulation 

The  finite  element  method,  in  conjunction  with  the  least  squ  es  method 
of  weighted  residuals,  is  used  herein  to  solve  numerically  the  small  perturba- 
tion transonic  equation.  In  this  approach,  a set  of  locally  defined  trial  func- 
tions with  undetermined  parameters  is  assumed  as  the  approximate  solution, 
and  the  integral  expression  for  the  square  of  errors  committed  by  the  approxi 
mate  solution  is  formulated.  Then  the  integral  of  square  errors  in  the  entire 
domain  is  minimized  with  respect  to  the  undetermined  parameters  to  yield  a 
system  of  algebraic  equations.  In  actual  computations,  the  minimization 
process  is  performed  at  element  level  and  then  an  assembling  process  is 
invoked  to  obtain  the  system  of  algebraic  equations.  It  is  these  processes 
which  enable  us  to  take  care  of  the  supersonic  pocket  with  convenience  and 
success.  This  special  assembly  technique  is  described  in  Section  3.3-. 

Written  in  the  usual  manner  with  repeated  indices  implying  summation, 
the  approximate  solution  has  the  following  form 
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0 = N.  0. 
i 1 


(9) 


in  which  N/ s represent  the  shape  functions  and  0 . ' s are  the  undetermined 
parameters  at  nodal  points.  The  residual  then  becomes 


R = j [ 1 - M2  - M2  (1  + y)  N,  0 ,1  N.  + N.  ,rir  j 0 . 
( L oo  oo  " k,x  r k J j ,xx  j.yy  Jrj 


(10) 


from  which  an  integral  expression  for  the  square  errors  is  obtained  as 


If1 
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Upon  the  minimization  of  X with  respect  to  the  undetermined  parameters,  a 
system  of  algebraic  equations  in  the  following  form  is  obtained 


S..  0 =0 


(12) 


Herein  the  banded  system  matrix  S. . is  defined  as 

ij 


S.. 
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. Q.  dA 


(13) 


with  and  CT  determined  by  the  following  expressions 

Q.  = [ 1 - M2  - M2  (1  +y ) N.  0.1  N.  + N. 

J l °o  oo  k,x  rk  j j,xx  j, 


yy 


(14) 


and 


P = Q.  - (1  + y)N.  0,  N. 

i i oo  ' k,xx  rk  1,3 


(15) 


As  stated  earlier,  the  system  matrix  is  obtained  by  combining  appropri 
ate  contributions  from  all  the  elements.  The  element  matrices,  in  turn,  are 
evaluated  effectively  by  numerical  integration  to  avoid  the  tedious  and  error 
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prone  algebraic  manipulations.  The  one  presently  used  is  a 7 point 
Gaussian  quadrature,  which  can  integrate  exactly  a quintic  polynomial  in 
the  element.  It  is  to  be  noted  that  although  the  system  of  equations,  Equation 
(12),  is  homogeneous;  it  does  possess  a non -trivial  solution  once  the  boundary 
condition  on  the  airfoil  is  imposed. 

3.2  Element  Description 

The  basic  element  used  in  the  present  program  is  the  nonconforming 
cubic  triangular  element  developed  by  Bazeley  et  al.  (Ref.  2).  Also  used  in 
the  program  are  quadrilateral  and  trapezoidal  elements  constructed  from 
these  triangular  elements.  These  three  types  of  elements  can  be  mixed  and 
used  freely  in  the  entire  flow  region  except  that  only  trapezoids  should  be 
used  to  cover  adequately  the  supersonic  region  in  order  to  enact  the  special 
treatment  discussed  later. 

The  basic  triangular  element  is  shown  in  Fig.  2,  which  at  each  vertex 
has  the  function  itself  and  its  two  first  derivatives  (velocity  components)  as 
undetermined  parameters.  This  type  of  element  was  adapted  mainly  because 
boundary  conditions  of  both  Dirichlet  and  Neumann  types  can  be  imposed  with 


Figure  2 - Triangular  Element  with  Undetermined  Parameters 
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equal  convenience.  In  addition,  because  velocities  at  nodes  are  treated  as 
primary  unknowns,  secondary  unknowns  such  as  local  Mach  number,  and 
pressure  coefficients,  etc.  can  be  computed  directly  without  resorting  to 
numerical  differentiation  and  thus  assuring  higher  accuracy.  Furthermore, 
the  use  of  higher  order  elements  can  usually  improve  computational  effi- 
ciency as  evidenced  in  most  finite  element  analyses. 

In  the  element,  the  approximate  solution  is  assumed  as 

= N.  (j).  (i  = 1 to  9) 

r l v l 

in  which  <j  's  are  the  nine  undetermined  parameters  of  <f)  and  its  first  deriva- 
i 

tives  at  the  nodal  points  as  shown  in  the  figure  and  N.'s  are  the  corresponding 
shape  functions,  which  are  expressed  in  terms  of  area  coordinates. 

Defined  in  the  following  are  these  shape  functions  and  their  first  and 
second  derivatives. 


Letting 


a.  = x.y  -■  x,  y. 
i J k k7j 

b.  = y.  - y, 

i j 7 k 

c.  = x.  - x. 
i k J 

A = area  of  triangle  1-2-3  = (b.c^  - b^c.)/2 

J J 

a = 0.5  (<=k-c.) 

P = 0.5  (b.-bk) 

Hx  = bi^k  + bj  Vi  + bk^i^  j 

Hy  = Vj^k  + Cj  Vi  + CkVj 

Hxx^^iVk^AVWj) 

Hyy  = 2<  ^ icjck  + ^ jCkCi  + ViCj* 


With  i = (l  .2,  3).  j = (2,  3,1),  k = (3.1,2),  then  one  has 
for  jf—  (1,4,7),  i=(l  ,2,3), 


Nf  = q (3-2£.)  + 2H 
Ni,x  = 2X  [6bi&*(1_Ci)  + 2Hx] 

N|.y=2X  [ i)  + 2Hy  i 

N<,xx  "TTTJ  K'1-2^'  * 2Hxx| 


l'*  - — 

l.Y 

2A  l 

T — 

1 

XX 

(2A)2 

J — 

1 

H - 

t.  vy 

(2A)2 

73  K<‘-2^  + 2Hyy] 


for  /-( 2,5,8),  i=(l,2,  3) 

r 2 , _ 


n<  = wnl*> + aH 

Ni,x  = 23l2bi^(CkS-CJ;k)  + 2A^+“Hx] 

Ni.y^KWci^)  + QHy] 


l,  XX 


TTK«ek‘j-cjtk>+4bi‘2A»it,<H»l 


J.  = — ^ [2c2(c,  - c.;,  ) + aH  1 

t.  yy  (2A)2  t i k j j k yyJ 


for  l=(3,  6,  9),  i=(l,  2,  3) 


■rfjt  »- 


VWk- Vj,  + PH 

N<.x=2i[2bi^i<bi^-bkS)  + PHx| 

Ni,y  = ^:K^(bj^-bkV  + 2ACi+pHyl 
N/.30c  = ^rK'Vk-bkV  + pHxxl 

N<.yy  * ^H<Vk  ‘ Vj>  + + PV! 


The  undetermined  parameters  are  arranged  in  the  order  of  ^ u , v , <*,  u 


v2>  0 3 , u^,  and  v^. 


Quadrilateral  and  trapezoidal  elements  as  shown  ii  Fig.  3 are  also  used 
in  the  present  program,  the  former  in  purely  subsonic  flow  region  to  save 
data  input  and  the  latter  in  the  mixed  and  supersonic  flow  region.  The  ele- 
ment matrix  for  a quadrilateral  is  obtained  by  combining  appropriately  the 
matrices  for  two  triangles,  while  matrix  for  a trapezoidal  element  is 


a.  Quadrilateral  Element 


b.  Trapezoidal  Element 


Figure  3 - Elements  Constructed  from  Basic  Triangles 


. E— . 


obtained  by  combining  and  averaging  contributions  from  four  triangles,  two 
left-running  and  two  right -running  triangles.  The  averaging  process  is  de- 
signed to  remove  the  bias  effects  inherent  in  the  quadrilaterals  presently 
used.  In  any  e'  ent,  the  present  program  is  coded  in  such  a way  that  the  com- 
puter is  given  the  jurisdiction  over  which  type  of  element  to  use,  according 
to  the  flow  behavior  in  the  element.  Of  course,  it  does  no  harm  to  use  trape- 
zoids where  the  usage  of  quadrilaterals  is  considered  adequate,  except  that 
computation  time  will  increase  slightly. 

3.3  Special  Considerations  for  Supersonic  Region 

As  is  well  known,  the  finite  element  method  is  a powerful  tool  for  solv- 
ing problems  governed  by  elliptic  equations.  However,  for  problems  governed 
by  either  parabolic  or  hyperbolic  equations,  the  conventional  finite  element 
assembly  technique  must  be  modified  so  that  the  proper  influence  of  solution 
in  time  (or  time-like  direction)  is  considered.  For  the  present  problem, 
the  x-axis  is  a time-like  direction  which  implies  that  the  solution  at  the  up- 
wind station  will  determine  the  solution  at  its  downwind  station  but  not  the 
contrary.  This  consideration  leads  to  the  well  known  upwind  influence  finite 
difference  operator  (backward  finite  difference).  An  equivalent  technique  in 
finite  element  analysis  has  been  devised  and  applied  to  solve  the  transonic 
flow  equation. 

Consider  the  rectangular  element  sketched  below  with  upwind  station  I 
and  downwind  station  II,  each  having  two  nodal  points.  With  the  element  type 
chosen,  the  element  matrices  can  be  constructed  in  the  usual  manner.  How- 
ever, before  assembling  the  element  matrices  into  the  system  matrix,  the 

coefficient  of  the 

1 

Flow  Direction 
— ►- 

2 

I I 

! I 

i I 

i II 
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first  term  in  the  transonic  equation  should  be  evaluated 


C - 1 - (1  + y)  u (16) 

The  sign  of  the  coefficient  C being  positive,  zero,  or  negative  will  define  the 
equation  as  elliptic,  parabolic  or  hyperbolic.  If  C is  non-positive  for  all 
nodes  in  the  element,  the  rows  representing  the  improper  downwind  influence 
on  solution  at  an  upwind  station  are  ignored  during  assembly.  This  feature 
is  taken  care  of  automatically  in  the  program,  requiring  only  a little  care 
with  the  nodal  arrangement  in  the  element.  In  the  anticipated  supersonic  re- 
gion, element  node  points  should  be  arranged  in  the  order  as  indicated  in  the 
above  sketch,  starting  with  the  upper  left  corner  node  and  proceeding  in  counter- 
clockwise direction.  On  the  other  hand,  if  the  sign  of  C is  positive  at  any  of 
the  four  nodes,  no  special  assembling  it  invoked.  As  stated  earlier,  only  tra- 
pezoidal elements  are  to  be  used  in  the  anticipated  supersonic  and  mixed  flow 
region,  while  triangular  elements  and  quadrilaterals  consisting  of  only  two 
triangles  are  considered  unsuitable  because  of  their  bias  nature.  However, 

these  two  latter  types  of  elements  can  be  used  effectively  in  the  subsonic  flow 
region. 

3.4  Boundary  Conditions  and'Iterative  Procedures 

As  stated  earlier,  the  imposition  of  boundary  conditions  for  the  present 
problem  can  be  carried  out  conveniently  because  the  elements  presently  used 
have  function  value  and  two  first  derivatives  as  primary  unknowns.  The  asso- 
ciated boundary  conditions  thus  can  be  treated  as  the  essential  type,  i.  e. , having 
prescribed  values.  Standard  finite  element  methodology  is  therefore  followed 
by  assembling  first  an  unconstrained  problem  and  then  modifying  the  matrix 
equations  accordingly.  More  specifically,  on  the  outer  boundary  the  flow  is 
assumed  to  be  undisturbed  and  hence  the  corresponding  unknown  parameters 
are  set  equal  to  zero.  On  the  line  of  symmetry,  the  velocity  component  v is 
also  set  equal  to  zero.  On  the  airfoil,  the  nonlinear  form  of  flow  tangency  on 
the  airfoil  is  imposed  by  replacing  the  algebraic  equation  for  v by 


in  which  u.  and  v^  are  unknown  parameters  at  node  "i"  on  the  airfoil.  When 
the  linear  form  is  desired,  the  corresponding  equation  becomes 


- & 
i dx 


(17a) 


and  is  applied  along  the  chordline.  The  selection  between  these  two  approaches 
is  made  through  a control  key  in  the  input. 

In  the  present  program,  no  special  efforts  have  been  devoted  to  treat  the 
leading  edge  and  trailing  edge  singularities.  Any  rigorous  approach  should 
consider  also  the  invalidity  of  small  perturbation  theory  in  these  regions, 
which  is  certainly  beyond  the  scope  of  the  present  study.  Nevertheless,  these 
regions  are  relatively  small  and  the  flow  is  usually  subsonic  in  nature,  which 
implies  any  error  committed  is  generally  localized  and  becomes  less  impor- 
tant elsewhere.  For  non-blunt  airfoil  such  as  a 6%  thick  circular  arc,  using 
half  of  the  actual  airfoil  slope  at  each  respective  place  in  the  computation  seems 
to  work  well.  However,  for  blunt  airfoil  such  ao  NACA  64  A006,  a finite  slope 
near  the  leading  edge,  say  1%  behind  the  airfoil  leading  edge,  could  be  used  in 
computing  the  vertical  velocity  component  at  that  point.  A still  better  approach 
is  to  use  very  fine  elements  in  that  region  and  treat  the  leading  edge  as  a singu- 
lar point.  In  practical  computations,  however,  an  element  with  its  area  ap- 
proaching zero  might  create  a new  numerical  problem  because  all  shape  func- 
tion derivatives  involve  element  area  as  divisor.  Therefore,  caution  must  be 
taken  in  using  extremely  small  elements. 

With  the  equations  properly  assembled  and  boundary  conditions  imposed, 
the  system  of  nonlinear  algebraic  equations  is  finally  solved  by  iterative  pro- 
cedures in  the  form 

S. . (6)  d(.n)  = L.  (18) 

ij  y rJ  i 
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in  which  0 is  a relaxation  factor  in  the  range  0 < 0 < 1.  For  subsonic  flows, 

0 is  usually  taken  as  unity;  for  flow  in  the  transonic  regime,  however,  it  has 
been  found  that  the  use  of  an  under  relaxation  factor  is  retired.  The  value  of 
0 should  decrease  from  unity  for  barely  critical  flow  to  0.4  or  even  0.3  for 
supercritical  flow.  The  optimum  underrr elaxation  factor  generally  depends 
on  the  freestream  Mach  number  and  the  mesh  being  used.  At  the  present,  no 
systematic  investigations  have  been  done  in  this  regard,  and  therefore  its  value 
should  be  selected  by  numerical  experimentations.  Generally  speaking,  a 
smaller  under-relaxation  factor  will  make  the  solution  more  stable  but  at  the 
same  time  tends  to  slow  down  the  rate  of  convergence. 

Equation  (18)  is  to  be  solved  subject  to  certain  prescribed  convergence 
criterion.  The  one  presently  used  is  that  the  relative  change  of  local  Mach 
number  between  two  consecutive  iterations  should  be  less  than  a prescribed 
small  number  for  all  the  nodes  in  the  flow  field,  that  is 


Numerical  experimentations  indicate  that  a solution  with  adequate  accuracy 
is  generally  obtainable  with  r in  the  range  0,01  <_  e <_  0.001. 


4.  PROGRAM  DESCRIPTIONS  AND  USAGE 


As  stated  earlier,  the  present  computer  program  is  capable  of  analy- 
zing steady  transonic  flow  over  nonlifting  thin  airfoils  by  solving  the  small 
disturbance  but  nonlinear  transonic  potential  equation.  The  numerical  al- 
gorithm is  based  on  utilizing  the  finite  element  technique  in  conjunction  with 
the  least  squares  method  of  weighted  residuals.  To  solve  a particular  tran- 
sonic flow  problem,  an  appropriate  finite  element  mesh  must  first  be  set  up, 
followed  by  using  the  numerical  procedures  summarized  in  Section  3.  Thus 
the  present  programs  are  separated  into  two  parts  — the  first  part  which  gen- 
erates the  necessary  mesh  information  from  a limited  number  of  input  cards 
and  the  second  part  which  carries  out  the  analysis  with  element  mesh  gener- 
ated in  the  first  part.  By  doing  so,  the  generated  mesh  can  be  fully  inspected 
for  its  correctness  prior  to  the  analysis,  and  storage  required  in  the  second 
program  can  be  accordingly  determined  to  suit  each  particular  problem.  These 
two  programs  are  designated  STRANL-I  and  STRANL-H  (Steady  Transonic 
Flow  by  Lockheed,  parts  I and  II,  respectively)  and  are  described  in  more  de- 
tail in  the  following. 

4.1  STRANL-I 

4.1.1  Scope  and  Flow  Chart 

This  program  reads  in  a limited  number  of  input  cards  and  generates  mesh 
information  such  as  element  nodes,  nodal  coordinates,  boundary  nodes  and  air- 
foil slope,  etc.,  to  be  used  in  the  second  program.  Additionally,  the  present 
program  has  two  options  — one  for  renumbering  the  nodes  to  yield  a smaller 
bandwidth  and  the  other  to  plot  the  generated  mesh  for  quick  inspection.  A 
schematic  flow  chart  of  the  program  is  shown  in  Fig.  4. 


Read 

TITLE 

IOPT  (Control  Keys) 

Call 

MESHBW 

(Generate  mesh  infor- 
mation from  a limited 
number  of  input  cards.) 

Read 

NIDS,  NID,  VAF 

(Boundary  nodes  and 
airfoil  slope) 

IOPT(II)  = I 


If 


NJNT  (I)  = I 

(Retain  original 
node  numbering) 


(Renumber  nodes 
to  reduce  bandwidth) 


IOFT(IZ)  = 1 


Print  Output  Data 


C/.ll  MPLOT 
(Plot  generated  mesh) 


Punch  Output  Data 
(or  store  on  tape) 


The  program  as  presently  dimensioned  requires  approximately  4Zg  won 
to  run  and  i an  aa  ommodale  the  following  maxima; 


• 200  elements 

• 180  nodal  points 

• 50  nodes  for  each  type  of  boundary  conditions 

• 10  nodes  connected  to  any  node  in  the  mesh 


4.1.2  Description  of  Variables 


TITLE 

IOPT 


NDEL 

X 


Y 

NJNT 


MEMJT 
NEW  JT 
JOINT 
JMEM 


NIDS 


NID 


VAF 


NEM 

NPM 


Array  used  to  describe  the  problem 
under  consideration. 


Array  containing  the  option  keys  which 
are  activated  when  read  in  as  "1."  Pre- 
sently there  are  two  keys  in  the  program, 
IOPT  (11)  for  node  renumbering  option 
and  IOPT  (12)  for  mesh  plot  option. 


Array  containing  element  node  points. 

Array  containing  the  x-coordinate  of 
nodal  points. 


Array  containing  the  y-coordinate  of 
nodal  points. 


Array  containing  the  new  nodal  numbers 
if  node  renumbering  option  is  executed; 
otherwise  it  retains  the  original  nodal 
numbering. 


Arrays  to  provide  adjustable  dimensions 

used  in  subroutine  OPTM  for  node  renumbering. 


Array  containing  the  actual  number  of 
boundary  nodes  for  far  field,  on  the  line 
of  symmetry,  and  along  the  airfoil  surface. 

Array  containing  the  nodal  numbers  for 
each  type  of  boundary  nodes  mentioned 
above. 


Array  containing  the  airfoil  slope  for  nod 
on  the  airfoil  surface. 


es 


Maximum  number  of  elements  allowed. 
Maximum  number  of  nodal  points  allowed. 


a!  i'LiifaJ.-.  - L . 
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NCN 


L 1 -1*1-  M*-  


NEMN 

NELS 

NPS 

MAXN 


NBW 


Maximum  number  of  nodes  connected  to 
any  node  in  the  mesh. 

Product  of  NPM  with  NCN. 

Total  actual  number  of  elements. 

Total  actual  number  of  nodal  points. 

The  maximum  difference  in  nodal  num- 
bers between  two  connected  nodes. 

Full  bandwidth  of  the  resulting  system 
matrix. 


4.1.3  Subroutines 

MESHBW:  This  subroutine  is  called  to  generate  mesh  information  in- 

cluding the  following:  array  of  element  nodes,  total  number  of  elements,  x- 
and  y-coordinates  of  nodal  points,  and  the  maximum  difference  in  nodal  num- 
bers between  two  connected  nodes.  In  generating  this  information,  a small 
amount  of  input  cards  is  required,  which  will  be  described  under  Section  4.1.4 


OPTM:  This  subroutine  is  called,  when  desired,  to  renumber  the  nodal 

points  so  as  to  yield  as  small  a bandwidth  as  possible.  The  algorithm  first 
establishes  a table  containing  the  nodal  numbers  connected  to  each  of  the  no- 
dal points,  and  then  considers  each  node  in  turn  as  the  origin  of  the  new  num 
bering  system  to  search  for  one  resulting  in  a smallest  bandwidth.  In  conse- 
quence this  subroutine  generates  an  array  of  pointers  relating  the  original 
numbering  scheme  to  the  new  numbering  scheme.  This  array  is  to  be  used 
in  constructing  the  matrix  equation,  and  is  used  again  to  refer  back  to  the 
original  numbering  at  output  time.  Using  this  method,  one  need  not  worry 
about  the  bandwidth  problem  and  is  unaware  that  any  renumbering  has  taken 

place. 


MPLOT:  This  routine  plots  the  mesh  generated  together  with  the  origin- 

al node  numbering,  using  the  CALCOMP  plotter.  Correctness  of  the  mesh  can 
thus  be  checked  quickly  prior  to  running  the  second  program. 
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4.1.4  Input  and  Output 


/ 


All  input  and  output  are  referred  to  the  original  numbering  scheme,  with 
input  in  the  form  of  punch  cards,  and  output  in  the  form  of  printout  and  punch 
cards.  If  subroutine  MPLOT  is  called,  the  mesh  is  also  plotted  for  the  con- 
venience of  inspection. 

Input 

Input  cards  to  this  program  should  be  prepared  and  provided  in  the  order 
described  below. 


P - 


: 


i 

i 


V 


A TITLE  CARD  (7A10,  A2) 

Col.  1-72  Description  of  the  problem  under  study 

B.  OPTIONS  CARD  (4012) 

Col.  22  Punch  "1"  if  node  renumbering  is  desired 

Col.  24  Punch  "1"  if  mesh  plot  is  desired 

C.  ELEMENT  CARDS  (1615) 

These  are  cards  supplied  to  generate  element  nodal 
numbers  in  groups.  One  card  is  needed  for  ea^h 
group  of  elements  whose  nodes  are  related  in  a 
regular  pattern.  The  number  of  elements  in  a group 
can  vary  from  one  to  any  positive  integer  number, 
depending  on  how  the  nodes  are  numbered  origin- 
ally. There  can  be  as  many  cards  as  required  to 
generate  the  element  nodes,  and  a blank  card  is 
added  at  the  end  to  terminate  the  process. 


Col. 

Description 

1-5 

Number  of  nodes  per  element.  Punch  3 for 
triangular  elements,  and  punch  4 for  quadri- 
lateral and  trapezoidal  elements. 

6-10 

Number  of  elements  in  one  direction  (called 
first  direction). 

11-15 

Number  of  elements  in  the  other  direction 
(called  second  direction). 

16-20 

Increment  of  nodal  numbers  in  the  first 
direction. 

21  -25 

Increment  of  nodal  numbers  in  the  second 
direction. 

26-30  First  nodal  number  in  the  element 

31  -35  Second  nodal  number  in  the  element 

36-40  Third  nodal  number  in  the  element 

41-45  Fourth  nodal  number  in  the  element 
("0"  for  triangles) 

As  stated  before,  element  nodes  are  presently  ordered  in  the  counterclock- 
wise direction,  and  the  upper  left  corner  node  should  be  taken  as  the  first 
nodal  point  in  an  element  located  in  the  anticipated  supersonic  region.  This 
is  required  to  enact  the  assembling  process  correctly  in  the  supersonic  pocket. 

For  example,  to  generate  element  nodes  for  a mesh  shown  below,  three 
input  cards  are  required:  one  card  for  the  first  12  elements,  second  card  for 
the  next  two  elements,  and  the  third  card  for  the  triangular  element.  The 
three  cards  are: 


4 4 3 4 1 2 1 5 6 

4 1 2 0 1 18  17  21  22 

3 1 1 0 0 20  19  23  0 

(Ended  by  a blank  card  when  all  elements  are  covered) 


4 8 


I®  I © 

© / © 


© © 


u 


r 


r 


i 


■ 


I 


I 


1 


1 

4 


D CARD  FOR  TOTAL  NO.  OF  NODES  (15) 

Col.  1-5  Total  number  of  nodes  for  the  entire 


flow  region. 

E.  NODE  COORDINATE  CARDS  (315,  5F10  0) 


The  x-  and  y-coordinates  of  nodes  are  also  generated 
in  connected  groups.  A connected  group  of  nodes  are 
those  falling  on  a straight  line  with  a constant  ratio  of 
nodal  distances  and  constant  increment  in  nodal  num- 
bers as  well.  The  intermediate  nodal  points  can  thus 
be  generated  by  linear  interpolation.  Again,  the  pro- 
cess is  terminated  by  encountering  a blank  card. 


Col. 

1 -5 
6-10 
11-15 


Description 

Node  number  of  the  beginning  nodal  point. 
Total  number  of  nodes  on  the  line. 


Increment  of  nodal  numbers  between  two 
consecutive  nodes. 


16-25 


Ratio  of  distances  between  three  conse- 
cutive nodes. 


26-35 
36-45 
46-55 
56  -6  5 


x-coordinate  of  the  oeginning  node, 
y-coordinate  of  the  beginning  node  . 
x-coordinate  of  the  end  node, 
y-coordinate  of  the  end  node. 


For  example,  to  generate  the  coordinates  for  the  nodes  depicted  below  with  a 
constant  distance  ratio  of  1.2,  the  input  card  should  read  as 


2 5 4 1.2  x2  y 2 *18  y18 


14 
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F.  CARD  FOR  NO.  OF  NODES  ON  BOUNDARIES  (1615) 


CoL 

1 -5 


De  sc  ription 

Number  of  nodes  on  the  outer  boundary 
(called  NFARF) 


6-10 

11-15 


Number  of  nodes  on  line  of  symmetry 
(called  NWAKE) 


Number  of  nodes  on  the  airfoil  surface 
(called  NBODY) 


G. 


H. 


CARDS  FOR  BOUNDARY  NODES  (1615) 

(a)  Nodes  on  the  outer  boundary  (NFARF  entries) 

(b)  Nodes  on  line  of  symmetry  (NWAKE  entries) 

(c)  Nodes  on  the  airfoil  surface  (NBODY  entries) 

CARDS  FOR  SLOPE  ON  AIRFOIL  (8F10.0) 


Airfoil  slope  for  nodes  on  the  airfoil  surface  (NBODY 
entries) 


Output 


The  output  from  this  program  is  in  the  form  of  printout,  punch  cards, 
and  also  a plot  if  the  mesh  plot  option  is  invoked. 


The  following  items  are  printed  in  the  order  as  described: 


Title  of  the  problem  under  study 


Control  keys  specified  in  the  options 
card 


Total  number  of  elements,  number 
of  nodal  points,  and  the  full  bandwidth 

Element  numbers  and  element  node 
points 


Nodal  numbers  and  their  correspond- 
ing x-  and  y-coordinates 


Nodes  on  each  type  of  boundaries 
Slope  for  nodes  on  the  airfoil  surface 


The  above  items  except  the  first  two,  together  with  the  NJNT  array  relating 
the  new  node  numbering  to  the  original  node  numbering,  are  then  punched 
and  saved  as  input  to  the  second  program  (STRANL-II). 
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4.2  STRANL-1I 


I 


1 

I 


l 

I 


4.2.1  Scope  and  Flow  Chart 

This  program  carries  out  the  analysis  following  the  numerical  procedures 
described  in  Section  3,  with  any  mesh  information  generated  and  supplied  through 
the  first  program.  In  summary  this  program  solves  the  small  perturbation 
transonic  potential  equation  via  the  least  squares  finite  element  technique.  The 
embedded  supersonic  pocket  is  taken  care  of  by  a proper  assembling  scheme 
and  the  resulting  system  of  nonlinear  algebraic  equations  is  solved  by  direct 
iterative  scheme.  The  present  program  has  three  options,  one  for  inputting 
nonzero  initial  guess  of  the  solution,  one  for  computing  in  the  same  run  solu- 
tions for  higher  Mach  number  (s)  using  results  computed  for  lower  Mach  number 
case(s),  and  the  last  one  for  selecting  the  form  of  boundary  condition  on  the 
airfoil.  A schematic  flow  chart  of  the  program  is  shown  in  Fig.  5. 

The  program  as  presently  dimensioned  requires  approximately  loOg  K 
words  for  the  program  itself  and  can  accommodate  the  following  maxima; 

• 200  elements 

• 180  nodal  points 

• 50  nodes  for  each  type  of  boundary  conditions 

• 84  in  full  bandwidth  of  the  resulting  matrix 
equations 

4.2.2  Description  of  Variables 

TITLE  Array  used  to  describe  the  problem  under 

consideration 

IOPT  Array  containing  the  option  keys  which  are 

activated  when  read  in  as  " 1."  Presently 
there  are  three  options:  IOPT(l)  for  con- 
tinuation for  high  Mach  number  cases  while 
using  existing  results  computed  for  lower 
Mach  number,  IOPT(2)  for  inputting  non-zero 
initial  guess  (referred  to  the  new  numbering 
system),  and  IOPT(3)  for  selecting  the  linear- 
ized boundary  condition  applied  along  the 
airfoil  chordline. 


Print  Computed  Results 


FPLOT 
on  the  Airfoil) 


Punch  S(I.  NEW) 


Converged  ? 


Update  Solution 


SLP,  RMLP 


w 11  , 


I 

t 


NOD 

S 

SLP 

X 

Y 

RML 

RMLP 

COF 

NJNT 

NIDS 


Array  containing  element  node  points 

The  resulting  system  matrix  with  main 
diagonal  terms  stored  in  the  column 
numbered  NHBW.  After  solving,  the 
column  numbered  NBW  contains  the 
solution. 

Array  containing  the  solution  obtained 
from  rexalation  procedures.  This  solu- 
tion is  used  in  generating  the  system 
matrix  to  solve  for  the  next  iterative 
solution. 

Array  containing  the  x -coordinate  of 
nodal  points. 

Array  containing  the  y-coordinate  of 
nodal  points. 

Array  containing  values  of  local  Mach 
number  at  nodal  points  for  the  current 
iteration. 

Array  containing  values  of  local  Mach 
number  at  nodal  points  for  the  previous 
iteration. 

2 2 

Array  containing  values  of  M^  + M^ 

(1  + y)u  evaluated  at  nodal  points. 

Array  relating  the  new  node  number- 
ing to  the  original  node  numbering 
schemes. 

Array  containing  the  actual  number  of 
boundary  nodes  for  far  field,  on  the  line 
of  symmetry,  and  along  the  airfoil,  re- 
spectively. 


NID 

VAF 

AR 

LR 

NEM 

NPM 


L, 


Array  containing  the  nodal  numbers  for 
each  type  of  boundary  nodes  mentioned 
above. 

Array  containing  the  airfoil  slope  for 
nodes  on  the  airfoil  surface. 

Array  to  store  values  for  use  in  sub- 
routine FPLOT. 

Another  array  to  be  used  in  the  sub- 
routine FPLOT. 

Maximum  number  of  elements  allowed. 
Maximum  number  of  nodal  points  allowed. 


} 

* 
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NCM 

NA 

NRM 

NELS 

NPS 

NBW 

NHBW 

NEQ 

ITGIV 

ZTEST 

RMAC 

SQMAC 

FI 

IRES 

POT 

UPT 

V 

U 

PRATIO 

CP 

DELM 


Maximum  full  bandwidth  allowed  in  the 
resulting  system  matrix. 

Maximum  number  of  nodes  allowed  on  the 
airfoil  surface. 

Maximum  number  of  equations  allowed. 

Total  actual  number  of  elements. 

Total  actual  number  of  nodal  points. 

Full  bandwidth  of  the  resulting  system 
matrix. 

Half  bandwidth  of  the  resulting  system 
matrix. 

Number  of  equations  in  the  resulting 
system  of  equations. 

Cut-off  number  of  iterations. 

Small  number  used  to  check  convergence 
based  on  the  relative  change  of  local  Mach 
number  between  two  consecutive  iterations. 

Freestream  Mach  number. 

Square  of  the  freestream  Mach  number. 

Under-relaxation  factor  with  0 < FI  < 1.0. 

Number  of  iterations. 

Perturbed  velocity  potential. 

Perturbed  velocity  component  in  the  x- 
direction. 

Perturbed  velocity  component  in  the  y- 
direction, 

Total  velocity  component  in  the  x-direction. 

The  ratio  of  local  static  pressure  to  stag- 
nation pressure. 

Pressure  coefficient. 

Change  of  local  Mach  number  between  two 
consecutive  iterations. 


4.2.3  Subroutines 

NEWK:  This  subroutine  is  called  in  the  main  program  to  generate  system 

matrix  for  the  entire  flow  field  by  assembling  appropriately  the  element 
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matrices.  The  elements  are,  in  general,  a combination  of  triangles,  quadri- 
laterals and  trapezoids.  For  a 4-node  element,  this  subroutine  also  deter- 
mines either  to  treat  it  as  a quadrilateral  or  as  a trapezoid,  depending  on  the 
local  behavior  of  the  governing  equation  being  elliptic  or  otherwise.  Subrou- 
tine EMQT  is  called  to  generate  the  element  matrices. 

EMQT:  This  subroutine  generates  element  matrix  for  elements  in  the 
form  of  triangles,  quadrilaterals  and  trapezoids.  Element  matrix  for  the 
later  two  types  is  obtained  by  combining  appropriately  (and  averaging  for  the 
last  type)  contributions  from  triangular  elements  which  are  in  turn  evaluated 
in  subroutine  EMTC. 

EMTC : This  subroutine  evaluates  the  element  matrix  for  a triangular 
element  by  numerical  integration.  The  one  presently  used  is  a 7 -point 
Gaussian  quadrature,  which  can  integrate  exactly  a quintic  polynomial  in  the 
element.  Shown  below  are  the  locations  of  these  Gaussian  points  with  their 
corresponding  weights  shown  on  page  31. 


0.05961587 

0.47014206 

0.79742699 

0.10128651 


DERV:  This  is  a subroutine  called  by  EMTC  to  evaluate  the  first  and 
second  derivatives  of  the  shape  functions  at  a Gaussian  point,  based  on  the 
expressions  listed  in  Section  3.2. 

BNDEQ:  This  is  an  equation  solver  for  banded  nonsymmetric  system  of 

algebraic  equations,  utilizing  the  direct  Gaussian  elimination  scheme.  In  this 
solver  the  system  matrix  is  arranged  in  a twisted  form  such  that  main  diagonal 
terms  are  stored  in  the  column  numbered  NHBW  and  the  right-hand  side  vector 
is  in  the  column  NBW.  After  solving,  the  column  numbered  NBW  contains  the 
solution  vector. 

FPLOT:  FPLOT  collects  and  stores  data  as  it  becomes  available,  and 

upon  signal,  produces  a printer  plot  in  practically  any  orientation  and  size. 

It  should  be  regarded  as  a general  purpose  output  routine  for  displaying  out- 
put data  in  graphical  form. 

Description  of  Parameters 

Size  of  the  main  storage  array,  and  it 
should  never  be  larger  than  800,  which 
-.or responds  to  400  points  to  be  plotted. 

Counter  initialized  — usually  IPNT  = 0 - 
in  the  calling  program.  . It  is  incremen- 
ted by  2 each  time  a new  data  point  is 
entered  in  AR. 


Ml 

IPNT 
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AR  Main  storage  array.  It  should  be  in  a 

dimension  statement  in  the  calling  pro- 
gram. For  example,  DIMENSION  AR(800). 

LR  An  array  of  bytes  used  to  hold  the  curve 

number.  It  should  be  dimensioned  for  Ml/2. 
The  type  declaration  LOGICAL*l  LR(400) 
should  be  in  the  calling  program. 

ISTOP  Flag  used  to  signal  that  all  data  has  been 

entered.  ISTOP  = 0 causes  data  to  be 
stored.  If  ISTOP  = -1  and  NC  = NCMAX 
the  program  immediately  branches  to  the 
plotting  section.  If  ISTOP  =1  and  NC'  = 
NCMAX  a data  point  is  stored  and  then  the 
plotting  section  is  entered. 

NC  Curve  number  for  which  data  is  being  en- 

tered. It  must  be  a positive  integer  less 
than  21. 

NCMAX  The  number  of  curves  to  appear  on  the  graph. 

This  is  the  largest  value  NC  will  have. 

VI  The  horizontal  coordinate  of  the  data  point  to 

be  plotted. 

V2  The  vertical  coordinate  of  the  data  point  to 

be  plotted. 


4.2.4  Input  and  Output 


Input 


The  bulk  of  input  to  this  program  is  provided  through  STRANL-I  in  the 
form  of  punch  cards.  To  run  several  cases  of  various  freestream  Mach  num- 
ber but  using  the  same  mesh  and  boundary  conditions,  three  additional  cards 
should  be  furnished  for  each  case  to  be  computed.  These  cards  should  be 
provided  in  the  following  order 

A.  TITLE  CARD  (7A10,  A2) 

Col.  1-72  Description  of  the  problem  under  study 

B.  OPTIONS  CARD  (4012) 

Col.  2 Punch  "l"  if  it  is  a continuation  in 

computation  from  one  Mach  number 
to  another. 


iur1.  . v., 


T” 

■ 'S| 
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Col.  4 Punch  " 1"  if  cards  for  nonzero  initial 

guess  is  furnished  in  the  input  data. 

If  so,  insert  them  following  those 
generated  from  STRANL-I. 

Col.  6 Punch  "1"  if  linearized  boundary  condi- 

tion along  the  chordline  is  desired. 

NOTE:  As  seen  in  the  flow  chart,  when  IOPT(l)  = 1,  the  last  two  options  become 
inoperative. 


C.  PROBLEM  PARAMETERS  CARD  (15,  3F10.0) 

Col.  1-5  Cut-off  number  of  iterations. 

Col.  6-15  Small  number  used  to  check  conver- 
gence based  on  the  relative  change  of 
local  Mach  number  between  two  con- 
secutive iterations,  usually  taken  in 
the  range  from  0.01  to  0.001. 

Col.  16-25  Freestream  Mach  number 


Col.  26-35  Under  relaxation  factor  with  the  value 
of  1.0  for  subsonic  cases,  and  decreas 
ing  for  supercritical  flows. 


For  the  first  case  of  the  run,  punch  cards  from  STRANL-I  are  placed  after 
the  above  cards.  For  subsequent  cases  using  the  same  mesh,  three  addition- 
al cards  as  described  above  are  required  for  each  case. 


Output 

The  output  from  this  program  include  printout,  referred  to  the  original 
numbering  system,  and  punch  cards  as  well  for  the  cases  whose  solutions 
are  convergent. 


The  printouts  are  in  the  following  order: 


• Title  of  the  problem  under  study 

• Convergent  limit 
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Control  keys  specified  in  the  options  card 

Total  number  of  element,  number  of  nodal  points,  and 
the  full  bandwidth  of  the  resulting  matrix  equations 

Element  numbers  and  element  node  points 

Nodal  numbers  and  their  corresponding  x-  and  y - 
coordinates 

Boundary  nodes  for  each  type  of  boundaries 

Slope  for  nodes  on  the  airfoil  surface 

For  each  iteration,  the  computed  results  at  all  nodal 
points  are  printed.  These  include  the  perturbed  velo- 
city potential,  total  velocity  components  in  the  x-  and 
y -directions,  the  value  of  M^  (1  + 2.4  u),  local  Mach 
number,  the  ratio  of  local  static  pressure  to  stagnation 
pressure  coefficient,  and  finally  the  change  of  local 
Mach  number  between  two  consecutive  iterations. 

For  each  iteration,  the  value  of  along  the  airfoil 
is  displayed  graphically. 


5.  SAMPLE  CASES 


Two  sample  cases  are  presented  herein  to  demonstrate  the  usage  of 
the  computer  programs  described  in  the  previous  section  regarding  input 
and  typical  output  from  the  programs.  The  two  problems  are  flow  over  a 
NACA  64  A006  airfoil  analyzed  with  a coarser  mesh  and  flow  over  a 6%  thick 
circular  arc  computed  with  a finer  mesh  where  the  nodes  are  renumbered 
to  reduce  the  bandwidth. 

In  computing  a transonic  flow  field  for  a given  airfoil  shape  with  the 
present  programs,  the  following  procedures  are  generally  followed: 

• A desired  mesh  is  sketched  with  each  node  assigned  a 
number. 

• Based  on  the  mesh  sketched,  an  appropriate  set  of  input 
cards  is  prepared  and  supplied  to  STRANL-I  to  generate 
the  required  mesh  information  in  the  form  of  punch  cards. 

• The  above  punch  cards,  with  three  additional  cards  for  each 
case  and  cards  for  non-zero  initial  guess  if  available,  are 
used  as  input  to  STRANL-II  to  compute  the  flow  field. 

• During  the  execution  of  STRANL-II,  results  for  each  iteration 
are  printed,  and  for  solution  satisfying  convergent  criterion 
the  solution  is  saved  in  the  form  of  punch  cards  for  possible 
later  use. 

The  following  paragraphs  describe  in  more  detail  the  input  and  output 
for  the  two  sample  problems. 

5.1  Flow  Over  a 6%  Thick  Circular  Arc  Airfoil 

This  problem  was  analyzed  using  the  mesh  shown  in  Sketch  1 which 
consists  of  154  elements  together  with  170  nodes.  Note  that  a relatively 
regular  mesh  is  adapted  in  the  anticipated  supersonic  region  while  the 
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remaining  elements  are  arranged  rather  arbitrarily.  Also,  the  nodal  points 
need  only  be  numbered  according  to  convenience,  as  the  node  renumbering 
option  will  be  invoked  to  renumber  them  to  yield  a smallest  bandwidth.  With 
the  present  mesh,  although  a little  more  effort  is  required  in  preparing  input 
cards  to  the  STRANL-I  program,  the  elements  are  considered  to  have  been 
used  more  effectively  than  a regular  mesh.  Summarised  below  are  the  input 
and  output  to  the  STRANL-I  and  STRANL-II  programs  for  the  present  problem. 


5.1.1  STRANL-I 
Input 

The  necessary  set  of  input  cards  is  listed  on  the  next  three  pages. 
Note  that  both  options  for  plotting  mesh  and  node  renumbering  are  to  be 
enacted,  and  the  final  mesh  will  have  nodes  on  the  actual  airfoil  surface 

rather  than  on  the  chordline. 
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Output 


As  stated  before,  output  of  this  program  includes  printouts,  punch  cards 
and  a plot  for  the  generated  mesh. 


Printouts  and  punch  cards  for  the  present  mesh  are  listed  on  the  next 
12  pages. 
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PRINT  OUTS  FROM  STRAIN-I  FOR  CIRCULAR  ARC  AIRFOIL 


STEADY  TRANSONIC  FLOW  — MESH  6 — 154  ELEMENTS  *170  NODES 


a-o-o-o 

-C-Q- 

C-Q” 

)-3  1 

-0-0 

40.  OF 

ELEMENTS3 

154 

NO. 

ELE.  NO 

. AND 

ELEMENT 

NODE 

1 

31 

21 

26 

0 

2 

6 

1 

2 

7 

3 

li 

6 

7 

12 

4 

16 

11 

12 

17 

5 

21 

16 

17 

22 

6 

26 

21 

22 

27 

7 

7 

2 

3 

8 

d 

It 

7 

8 

i i 

9 

17 

12 

13 

18 

Id 

22 

17 

18 

23 

11 

27 

22 

23 

28 

12 

8 

3 

4 

9 

1? 

13 

8 

9 

14 

14 

13 

13 

14 

19 

15 

23 

18 

19 

24 

16 

28 

23 

24 

29 

17 

9 

4 

5 

10 

18 

14 

9 

10 

15 

19 

19 

14 

15 

20 

20 

2h 

19 

20 

25 

21 

29 

24 

25 

30 

22 

26 

27 

32 

0 

23 

27 

28 

33 

32 

24 

28 

29 

34 

33 

25 

29 

30 

35 

34 

26 

32 

33 

36 

0 

FULL  BANDWIDTH*  34 
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NOOE 

NEW  NOOE 

X(I) 

Y4I) 

1 

167 

0. 

0. 

2 

168 

. 51 55  IE  +00 

0. 

3 

169 

. 9845  3E  + 00 

0. 

4 

170 

• 14115E+01 

0. 

5 

166 

• 180  0 0E  + 01 

0. 

6 

161 

0. 

.3dOOQE+Od 

7 

162 

• 51 83  8E  + 00 

•24272E+00 

a 

16  1 

. 390 luE+00 

.19060E+00 

9 

1^-4 

• 1419+E+Ul 

• 143L6E+L0 

10 

165 

• 18100E  + 01 

•lODOOE+OO 

11 

153 

0. 

. 700  OOE  + GO 

12 

154 

• 52124E+00 

.57399E+C0 

13 

155 

. 9955  7E  + 00 

.45931E  + U 

14 

156 

* 14272E+01 

* 35496E  + 00 

15 

157 

. 1320  3E  + 01 

.26300E+00 

16 

144 

0. 

.11000E+01 

17 

145 

. 5661 + E + 00 

. 896C7E+GJ 

Id 

146 

» 10  55  7E+01 

•71865E+00 

19 

147 

• 14  87  2E  + 01 

. 56429E+00 

2d 

148 

. 1860  li  E + 01 

• 430  0 0E  + 00 

21 

134 

0. 

. 15000E+GI 

22 

135 

•617+3E+00 

• 1223  3E  + 0 1 

23 

136 

.11324E+01 

. 100  46E  + 01 

24 

137 

• 1561 3E  + 01 

• 81670E+0J 

25 

138 

* 192Q3E+01 

• 660  0 0£  +03 

26 

122 

• 5003 JE+OG 

* 2030  0E  + C1 

27 

123 

• 3353  2E+00 

•16735E+G1 

28 

124 

* 13315E  +01 

. 13764E+01 

29 

125 

. 1592 JE  + 01 

. 110  50  E + LI 

30 

126 

• 20  200E  + 01 

. B6000E  + 00 

31 

133 

0. 

. 20000E+G1 

32 

111 

• 12000E+01 

. 200D0E+L1 

33 

112 

* 1535  0E  + 01 

• 16640E+C1 

34 

113 

• 1841 8E+01 

. 13582E+C1 

35 

114 

. 21200E  + 01 

. 10800E+01 

36 

ldl 

• 1300  0E  + 01 

• 200  C0E  + G1 

37 

102 

« 2u  40  6E  + 01 

• 15455E+01 

38 

103 

•22500E+01 

. 115Q0E+01 

39 

91 

. 2200UE+01 

• 200  0 OE  + Ol 

40 

92 

•22953E+01 

< 15329E+C1 

AOMl.  ,.aM 
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2380QE+Q1 
25  OQQE  + Ol 
25QQ0E+Q1 
25dO JE+Q1 
25000E+01 
27QJ7E+Q1 
2620 JE+Ol 
J2Q0UE+Q1 
29594E+01 
27500E+Q1 
38000E+Q1 
34640E+Q1 
J1582E+01 
2dd00£+01 
45000E+01 
43  64  7E  + Q1 
36635  E *01 
» 33  08  3E  + 01 
, 2980  0E+01 
» 50  0 00  E + 01 
. 50  000E+01 
.43826E+01 

• 38676E+01 
•34392E+Q1 

• 30  8Q3E+Q1 
. 50000E  + 01 
. -+4278E  + 01 
. 39299E+01 

• 34958E+01 
.3120GE+01 
. 5j  00  0 E + 01 

• 4478  2E  *01 
•39890E+01 
.35 493E+01 

• 31 50  0 E 4-21 
. 5G003E  + 01 
• 44845E+31 
•4L154E+01 
• 35885E+01 
• 32000E+01 

• 51  0 3 3 E + Ql 
. 44845E+01 
•40 154E+01 


.12200E+01 

• 20Q00E<  31 
.160732+31 
•12500E+01 

• 200002  + 01 
•15529E+01 
•1220024QL 
• 20000E+01 
•15455E+01 
•1150QE+01 
•200002+01 
•16567E+01 
.13443E+Q1 
.106002401 
•20000E401 
. 16735E4  01 
•13764E+01 

• 11060E  4Q1 
.860002400 
•20000E4C1 
.150002401 
.12299E401 
.100462401 
.81670E+03 
.66000E4CO 
.110002401 
.905202400 
.735722400 
.588282403 
.460002403 
.700002+00 
.573992+00 
.459312+03 
.354962+00 
.26000E+00 

• 303  OOE  + QO 
.242722+00 
.190602+00 
.143162+03 
.1Q0U0E+Q3 

0. 

0. 

0. 


84 

1 

85 

6 

06 

160 

87 

159 

08 

158 

09 

152 

90 

151 

91 

150 

92 

149 

93 

143 

94 

142 

95 

141 

96 

140 

97 

139 

98 

132 

99 

131 

100 

130 

101 

129 

102 

128 

103 

127 

104 

121 

105 

120 

106 

119 

107 

118 

108 

117 

109 

116 

110 

115 

111 

110 

112 

109 

113 

108 

114 

107 

115 

106 

116 

105 

117 

144 

118 

100 

119 

99 

120 

98 

121 

97 

122 

96 

123 

95 

124 

94 

125 

90 

126 

89 

. 35885E+01 
. 3£0  OkJE^Ol 
.2DOOQE  + 01 
. 2UOO0E+O1 

• 2aQ0QE+01 

• 20  500E+01 

• 20  500E+Q1 
.205O3E+O1 
. 2U53  OE  + Ol 
.21250E+O1 

• 2125  OE  + Ol 
.21250E+O1 
.2125 0E  + O1 
. 2125  OE  +01 

• 22  0 0 0 E + 01 

• 220  0 0E  + 01 

• 220  00  2+01 
. 220  0 0E  + 01 

• 22  0 DO  E + 01 

• 22  000E  + 01 
. 2275  0E+01 

• 22  75  QE+01 
.22  75  JE  + 01 
•22750E+01 

• 22  75 0E+O1 

• 22  75QE+01 

• 22  750E+O1 
. 25  50  OE  +31 

• 2350  OE  +01 
. 23  5 0 QE  +01 
. 23  50  J E + 01 
.23500E+01 

• 2 J500E  + 01 

• 2350  OE  + 01 

• 2425  OE  + Ol 
. 2425  0E  + O1 
. 2425UE+01 

• 2425  OE  + Ol 

• 2425  OE  + Ol 

• 2425  OE  + Ol 

• 2425  OE  + Ol 

• 25  00  3E  +01 

• 250  0QE  + Q1 


0. 

0. 

0. 

.1O0OQE+GO 
. 220  0 OE  + 00 
. 57200E-02 

• 100  30E  + 00 
.213802+00 
. 350D0E+C3 
* 131502-01 
• 11316E+00 
.23317* +03 

• 377 18E  + 03 
.55000* +03 
» 19220E-01 
.11742E+00 

• 235252  + 03 
.376682+00 
.546372+00 
.750002+00 

• 2394  02-01 
. 11418E  + 03 
• 22246E+C3 

• 35241E  +0  0 
.508342+00 
.695462+00 
. 920  0 OE  + 00 
. 27310E-GI 
• 12124E+QQ 
.233952+00 
•36921E+00 
• 53151E+00 
• 7262BE+00 
•960002+00 

• 29330  2 -Cl 
.127082+00 
.244382+00 
. 38515E+03 
.554062+00 

• 756762  + 03 
. 100  OOE  + Ol 
.30000E-01 
. 1297IE+  00 
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127 

88 

• 250  0 QE  +01 

• 24934c  + 08 

128 

87 

• 250  0 0E  + 01 

• 39290E+00 

129 

86 

• 25  4 0 0 E +0 1 

•555L8E+0Q 

130 

85 

• 25b  0QE+01 

• 77192E  + 00 

131 

84 

•25QQ3E+Q1 

• 102  0 QE  + 01 

132 

80 

• 25  75  OE  + Ol 

•2933QE-G1 

133 

79 

•25750E+01 

•127Q8E+QJ 

134 

78 

. 25  75  0 £ + 01 

• 24438E  +00 

133 

77 

• 25  75  QE  + Q1 

. 38515E+00 

1.36 

76 

•25750E+01 

. 554Q6E+00 

137 

75 

• 2575  0E  + 01 

•75676E+00 

138 

74 

Er750E+01 

• 10000E+01 

139 

70 

• c‘«5  0 GE  +01 

•27310E-01 

140 

69 

.2j5Q0£+Q1 

• 12124E+00 

141 

68 

• 2b  50  0E  + 01 

•23395E+00 

142 

67 

•2650JE+01 

• 36921E+G0 

143 

66 

•26500E+01 

. 53151E  + 00 

144 

65 

•26503E+01 

• 7262  8E  + 00 

145 

64 

• 2 650  0E  + 01 

• 960  OOE  + OO 

146 

60 

•27250E+01 

• 23940E-01 

147 

59 

• 2725  0 E + 01 

• 11418E  + 03 

148 

58 

•27253E+01 

• 22246E+03 

149 

57 

• 2725  0E+01 

• 35241E+00 

150 

56 

• 27  25  GE+01 

•50334E+00 

151 

55 

•27250E+01 

• 69546E+C0 

152 

54 

.27 250E+01 

. 92000E+00 

153 

49 

. 2300GE+01 

•1922QE-G1 

154 

48 

• 280  0 0 E+01 

• 11742E+03 

155 

47 

• 2800  3E+Q1 

•23526E+G0 

156 

46 

• 28000E+01 

• 37668E  + 00 

157 

45 

• 280  0 JE+01 

• 54637E  + 00 

158 

44 

• 23  30 3 E + Ql 

•750 OOE+OO 

159 

37 

• 2875  OE+Ol 

• 13150 E- 01 

163 

36 

• 2875  OE  + Ol 

• 11316E+00 

161 

35 

.2875  OE+Ol 

•23317E+00 

162 

34 

• 23  75  OE  + Ol 

• 37718E+10 

163 

33 

•23750E+01 

• 55Q0QE  + 0J 

164 

27 

• 29500E  + 01 

• 5720  QE-02 

165 

26 

•29500E+01 

• 100  30  E + 00 

166 

25 

• 29500E  +01 

• 21380E+00 

167 

24 

•29500E+01 

. 35000E  + 00 

168 

15 

• 30  000E  + 01 

0. 

169 

14 

.3200 JE+01 

• 10000E+U 3 

170 

13 

• 30  OOOE  + Ol 

•22000E+00 

Sfetrift* »•»  ?«■  m.1-  -» 
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5.1.2  STRANL-II  Program 
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With  the  present  mesh,  three  cases  of  freestream  Mach  numbers  are 
to  be  computed.  These  are 

M = 0.806  (Subcritical) 

00 

M = 0.861  (Barely  Critical) 

oo 

M = 0.909  (Supercritical) 

00 


Input 

Listed  below  are  typical  input  cards  to  the  STRANL-II  program  to 
compute  flow  field  for  the  above  cases. 


STEADY  STATE  SOLUTION  — CIWCULAW  AWC 
0 0 0 

b O.OOb  0*80 6 1.000 

( Insert  cards  generated  in  STRANL-I) 
STEADY  STATE  SOLUTION  — CIOCULAW  AKC 
1 

8 O.Ol  0*861  1.000 

STEADY  STATE  SOLUTION  — ClWCUtAW  AWC 
1 

12  0*01  0*909  0*500 


170  NODES  — M=0*506 


170  NODES  — M=0*b61 
170  NODES  — M=0*909 


The  above  cards  indicate  there  is  no  non -zero  initial  guess  furnished  in  order 
to  compute  the  case  M = 0.806,  while  cases  with  higher  Mach  nupnber  will 
use  results  computed  for  lower  Mach  number  as  initial  guess.  In  all  cases, 
the  nonlinear  boundary  condition  on  the  actual  airfoil  surface  is  to  be  used, 
as  suggested  by  IOPT  (3)  = 0 in  the  options  card. 


Output 

The  output  from  this  program  includes  printouts  and  punch  cards  for 
the  convergent  solution(s). 
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The  printouts  are  in  the  following  order: 


• Title  of  the  problem  under  study 

• Convergent  limit 

• Control  keys  specified  in  the  options  card 

• Total  number  of  elements  number  of  ^ points, 
and  the  full  bandwidth  of  the  resulting  matrix 
equations 

• Element  numbers  and  element  node  points 

• Nodal  numbers  and  their  corresponding  x-  and 
y -coordinates 

• Boundary  nodes  for  each  type  of  boundary 

• Slope  for  nodes  on  the  airfoil  surface 

• For  each  iteration,  the  computed  results  at  all  nodal 
r»r»i nt s are  printed,  together  with  a printer  plot  dis 
playing  the  pressure  coefficient  on  the  airfoil  surface. 

The  above  items  except  problem  parameters  and  computed  results 
have  also  been  printed  in  the  STRANL-1  program,  but  are  repeated  here  or 
completeness.  For  brevity,  only  result,  for  a typical  iteration  are  listed 

the  next  five  pages  for  reference. 

5.2  Flow  Over  a NACA  64  A006  Airfoil 

This  problem  was  analysed  using  the  mesh  shown  in  Sketch  2,  which 
consist,  of  120  quadrilaterals  with  150  nodes.  For  this  particular  mesh,  it 
is  apparent  that  nodal  numbering  along  the  shorter  (verticall  dtrectmn  yields 
the  smallest  bandwidth  and,  therefore,  no  nodal  renumbering  is  required. 

Since  the  present  mesh  has  a rather  regular  pattern,  less  input  car  s are 
required.  The  procedures  described  for  the  circular  arc  atr  ot  pro 
also  applied  here.  For  brevity,  only  the  input  cards  to  both  programs  and 
some  typical  printouts  from  STRANL-1I  are  listed  here  for  reference.  Again, 
the  nonlinear  boundary  condition  on  the  actual  airfoil  surface  ts  to  be  used  tn 
computations  for  all  the  cases. 
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Figure  7-  Freehand  Sketch  of  Finite  Element  Mesh 
for  NACA  64  A006  Airfoil 
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INPUT  CARDS  TO  STRANL-II  FOR  NACA  64  A006 


STt-ADY  TRANSONIC  FL0W--NACA64 A006 — MESH  5 — M=0«B2b 

(blank  card) 
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(Insert  cards  generated  in  STRANL-I) 

STEADY  TRANSONIC  FLOW — NACAb4A006--MESH  s — M=0.B75 
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STt-ADY  transonic  fLOW--naca64A006 — mesh  b — m=o.^oo 
i 

12  o »0 1 0.90  0.500 

STEADY  TRANSONIC  FLOW — NACA64A00& — MESH  b — M=0.OA0 
1 

16  0.01  Q. 940  0.400 
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Appendix  A 

FORTRAN  LISTING  OF  STRANL-I 
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Appendix  A 


A FORTRAN  listing  of  the  source  deck  for  the  STRANL-I  program  is 
presented  in  the  following  pages.  The  program,  as  presently  configured, 
required  42gK  words  to  execute. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

1 o 
'->0 


c 


120 


c 

c 

c 


1 JO 
1 40 
1 bo 


c 

c 

c 


dd.0 


2J0 


' y.  • 


PwOGRAM  MA  1 N(  i NPUT  = 6b  • OUTPUT  = 1 J 1 ,PUNCH=8b»  TAPEb=  1 NPUT  » 

i tami.6=ouTPu  r ) 

GLNLRATL  2~U  MLSH  CONSISTING  UF  TRIANGLLS  ANU  UUADR  1 LATLWALS 
MWt.Jt.lMr  PivJOKAM  HAS  TWU  OPTIONS 
1 Up T ( 1 1 ) = 1 RENUMBER  NODES  TU  YIELD  JMALLER  BANDWIDTH 
1 op  r ( 1 2 ) = 1 CALL  MPLUT  To  PLOT  oENERATEU  MESH 

THE  PROGRAM  A J PRESENTLY  O 1 MENS 1 ONLO  ALLOWS  THE  FOLLOWING  MAXIMA 
200  ELEMENTS *180  NODES. 10  CONNECTING  NODES  TO  ANY  NODE. 

SO  NODES  FOR  EACH  TYPE  OF  BOUNDARY  CONDITIONS 

dlveloplo  and  coded  by  stevens  chan  of  lockheed-huntsville. ai.a. 

except  SUBROUTINE  mplot  which  Is  coded  AND  supplied  by 

CART.  GERRY  VANNEuREN  OF  AFFDL.  WRIGHT  — PATT LRoON  AIR  FORCE  BASE 
DIMENSION  T I T L E { 8 ) * I OPT (20)  .NOEL (200.4)  * MEM JT ( 1800) 

DIMENSION  X ( i 80 ) > Y ( 1 80 ) . NJNT ( 1 80 ) .NEW JT ( 1 80 ) . JO  I NT ( 180). JMEm  < 1 80 ) 

DIMENSION  N1DS( 3) *N1D(3. 50) .VAF(bO) 

Lou  I V ALENCE  ( N I DS ( 1 ) . NF  ARE  ) . (NlDS(2)  .NWAk-E)  » ( N IDS  ( 3 ) . NBODY  ) 

DATA  NEM . NPM » nCN  » NEMN/200 . 180*  10. 1800/ 

READ  ANU  GENERATE  MESH  INFORMATION  .ETC. 

READ (5.805) < T I T LE < I ) . 1 = 1 . 8 ) 

1 F ( EOF ( 5 ) ) 2oOQ .50 
continue 

READ ( 5.820)  ( I OP  T l 1 ) . 1 = 1 . 20 ) 

CALL  MEohB w ( (MEM  * NPM  ♦ NDEL  . X . Y . IMEL J . NPw  . MAXN  ) 

READ  BOUNDARY  NUDES  AND  AIRFOIL  SLOPE 
READ  825.  NlDS 
DO  120  1=1.3 
NS  = N I DS  ( 1 ) 

READ  825.  ( NIDI  I t J ) . J = 1 .NS ) 

READ  830.  < VAF ( I ) « I = 1 .NBODY ) 

execute  node  renumbering  and  mesh  plot  options 

IF  ( I OPT (11)  .EU.  1)  GO  TO  140 
DO  130  1=1, nPS 
NJnT(  1 ) = 1 
GO  TO  150 

CALL  OPTM ( NPM  » NtM.NPS.NELS, 4 . NCN , NEMN , MEM JT . NEW JT .JOINT. 

1 JMEM, NDEL .NJNT .MAXN) 

CONT INUE 

IF  ( 10PT ( 1 2 ) . EO.  1 ) CALL  MPLOT ( X , Y » NDEL . NPM , NEM , NELS . NPS ) 

compute  full  bandwidth  and  print  output  data 

NBW=6* (MAXN+1 ) 

PRINT  DIO.  < T I TLE( I ) . 1= 1 ,8) 

PRINT  820.  ( IUPT(  1 ),  1=1,20) 

PRINT  920,  NELS.NPS.NBW 
PRINT  9 JO 
DU  220  N= 1 ,nELS 

PRINT  825,  N, ( NDEL ( N» J ) . J= 1 , 4) 

PRINT  9Jb 
DO  230  1 = 1 ,NPS 

PRINT  940,  1 »imJNT  ( 1 ) «X(  1 ) ,Y(  1 ) 

PRINT  951,  (NIDI  1*  1 ).  1 = 1 .NFARF) 

71 


T 


PM  I (SIT 
PM  I NT 
PMINT 


952. 

953. 
955. 


(N1DI2*I)«1=1 .NwAKE) 
( N 1 D < 3. I ) . 1 = 1 .NBODY) 
( VAF ( I ) . 1 = 1 « NBODY ) 


C 

c 

c 


PUNCH  OUTPUT  UA 1 A TO  BE  U8ED  AS  INPUT  TO  STMANL — II 


PUNCH 
PUNCH 
PUNCH 
DO  350 


025. 
025. 
040. 
I = l 


350 


(NIDI  I . J)  . J= I .NS) 

( VAF  < I ) . I = I .NBODY) 
( NJNT ( 1 ) . I = 1 . NPS ) 


005 
020 
825 
830 
040 
91  0 
920 


930 

935 

940 

951 

952 

953 
955 

2000 


NELS  » NPS  »NBW  « <NIDS(  I > . 1=  I .3) 

( I WDEL ( 1 .J)  . J= I .4)  . 1 = 1 .NELS) 

(XU)»Y(I),I  = I .NPS) 

.3 

NS=NIDS( I ) 

PUNCH  825. 

PUNCH  040. 

PUNCH  025. 

GO  I O I O 
FOMM AT  < 7A 1 0 . A2 ) 

FORMAT! 4012) 

FORMAT  (1615) 

FORMAT  ( OF  1 0 • 0 ) 

FORMAT  ( 1 POE  10 • 3 ) 

FORMAT! 1HI .2X • 7A10.A2// ) 

FORMAT  (I  HO, -NO.  OF  ELEMENTS='  • 14.  • NO. 

• FULL  BANDW1DTH=* . 14//) 

( 1 OELE • NO.  AND  ELEMENT  NODES'/) 

‘•OOLD  NODE'.'  NEW  NODE'.  6X . ' X ( I > ' . I 2X . • Y ( I ) • / ) 
1 16.  I 10 . 2E I 5.5 ) 

('ONOUES  AT  FARF 1 ELD • / ( 20  I 5 ) > 

UN  LINE  OF  SYMMETRY' /( 2015) ) 

ON  THE  A 1 MF  QlL'/(20l5>  ) 

ALONG  NODES  ON  A I MFO I L • / ( 8E I 5 . 5 > ) 


OF  NODES = • . 14 . 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

STOP 

END 


< • ONODLS 
( • ONoDtS 
( 'OSLOPE 


C 

c 

c 

c 


c 

c 

c 


100 


lEo 


124 


126 

129 


C 

C 

C 


130 


134 

135 


SUBROUTINE  mESHBWI NEM  «nPM«  NDEL » X,  Y «NELS « NPS » MAXN ) 


GENERATE  MEsH  information  including  element  node  array  * 
nodal  COORDINATES*  And  max.  difference  in  nodes 


1020 
1 150 


DIMENSION  NdEL{NEM,4> .X(NPM)  .Y(NPM) 
DIMENSION  I DUMP ( 4 ) »DUMP( 4 ) 

FORMAT  < 1615) 

FORMAT  (3I5.SF10.0) 


GENERATE  ELEMtNT  NODE  NUMBERS  AND  THEIR  MAXIMUM  DIFFERENCE 


DO  100  N= 1 .NEM 
DO  100  U= 1 .4 
NDE L < N « J ) =0 
NElS  = 0 
MAXN=0 

READ  10E0.  NPL.KEI .KEJ.KDI .KDJ. ( IDUMP( I ) , 1=1 ,4) 

IF  (NPE  . iU.  Q)  GO  TO  129 

Jl =-KDJ 

DO  124  J=1 «KEJ 

J1  = Jl  + KDJ 

J2  = Jl  - KDl 

DO  124  1 = 1. KEi 

NELS  = NELS  + 1 

J2  = J2  + KDl 

DO  124  K=1 .nPE 

NDEL(NELS.K)  = I DUMP ( K ) + J2 
DO  126  1=1, NPE 
DO  126  J=I ,NPE 

NDiFF=NDEL<nELS. II -NDEL < NELS, J) 

IF  (IABS(NDIFF)  ,GT.  MAXN)  MAXN= I Abb ( ND I FF > 

GO  TO  120 

IF  (NELS  «Gt»  NtM)  GO  TO  300 


GENERATE  NOdAl.  COORDINATES 


NPS 


0)  GO  TO  140 


3)  GO  TO  135 


READ  1 020 , 

READ  1150,  ND  * NN I , ND I , RATIO, ( DUMP ( I ) , I = 1 , 4 ) 
IF  (ND  «EQ< 

TEMP2  = 1.0 
IF  (NN1  .LT. 

Jl  = NNI  - 1 
TEMPI  = 1.0 
TEMP 2 = 0.0 
DO  134  1=1 , J 1 
TEMP 2 = TEMp 2 
TEMPI  = 

TEMPI  = 

TEMP2  = 

Jl  = ND 
TEMP3  = 0*0 
TEMP4  =0*0 
DO  137  1=1, NNI 
X( Jl  ) = DUMp( 1 ) 

Y< Jl > = DUMp( 2) 

Jl  = Jl  + NO I 
TEMP 3 = TiMp3  + 


+ TEMPI 
RATI 0*TfcMP 1 

(DUMP! 3) -DUMP ( 1 ) >/TEMP £ 
( DUMP l 4 ) -DUMP  < 2 ) ) /TEMP2 


TEMP3 

TEMP4 


TEMP  1 


73 


mu 


h'-'-ii.imrt'T-ii  ■ *r  I -•  . , 


TLMP4  = TEMp4  + TEMP2 
T1MP1  « KA  1 1 0*  T LMPl 
i w n.MPi:  * rat iu*TLMP£ 

OO  TO  I JO 

1A0  IF  < NPS  • OT • NPM)  GO  TO  400 

prInt^hhor  message  if  dimensions  for  elements  or  nodes  exceeded 

300  print  310*  nELS.NEM  .kbtm-i.Iai 

310  FORMAT  COTOO  MANY  ELEMENT  b*.bX.'  NELS=  • • I 4 • 5X  • • NEM-  • . I 4 ) 

CALL  exit 

400  PRINT  4 10*  NP  S » iMPM 

410  FORMAT  I'OToO  MANY  NODES'*  5X*  *NPS=  • * 14  *£>X*  NPM-  *14) 

CALL  EXIT 
END 


I 

* 


Wr  r-*?' 


L 

C 

L 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


SUBROUTINE  OPTMf  NPM.NEM.NPS.NELS.NPE .NCN.NEMN.MEMJT .NEW JT • 

1 JOINT. JMLM.NOEL.NJNT .MAXN) 

RENUMBER  the  NODES  TO  REDUCE  BANDWIDTH.  WITH  NJNT  CONTAINING  THE 
NEW  nodal  NjMDEKS  (LOCATION  IN  ARRAY  RERRLSENTS  THE  OLD  NODAL  NO' 
NRM-MAX*  no.  UF  MODES.  NPS=ACTUAL  no.  of  nodes 
NEM=MAX.  NO.  OF  ELEMENTS.  NELS*ACTUAL  NO.  OF  ELEMENTS 
NHe=MAX«  no.  OF  NODES  AMONG  ALL  THE  ELEMENTS 

ncn=est 1 mated  max.  no.  of  nodes  connected  to  any  node 
NDEL=ELEMENT  node  points 

MAxNsMAX.  DlFF.  BETWEEN  NODAL  POINTS  IN  ANY  ELEMENT 
MEMJT.NEWJT. JOINT. JMEM  ARE  INTERMEDIATE  WORKING  SPACE 

DIMENSION  NdEUNEM.NPE)  .NJNT(NPM)  , MEMJT  ( NEMN) 

DIMENSION  NEWJKNPM)  .JOINT  ( NPM  > , jMEM(NPM) 

ESTABLISH  NoDtS  connected  to  each  of  the  nodal  point 
DO  10  J=I.NPS 

io  jmlm(j>=o 

DO  60  J=1.N£LJ 
DO  50  1=1 «NPE 
JNT 1 aNDEL ( J . I ) 

IF  (JNT I .EQ.  0)  GO  TO  60 
JSuB=( JNTI-I ) *NCN 
DO  40  I I - 1 « NPE 
IF  (II  .EQ.  I)  GO  TO  40 
JJT  = NDEL( J« I I ) 

IF  ( JJT  .EG.  0 ) GO  TO  50 
MEM  I x JMEM ( JnT I ) 

IF  (MEMI  .EQ.  0)  GO  TO  30 
DO  20  II  I = 1 ,MEM1 
JP= JSU6+I I I 

IF  (MLMJT(JP)  .tQ.  JJT)  GO  TO  40 
20  CONTINUE 

30  JME.M(  JNT  I >=JMEM(  JNT1  ) + I 
NPoS  = JSUB+JmEM( JNT  1 ) 

MEMJT(NPOS)=JJT 
40  CONTINUE 
bo  continue 
60  CONTINUE 

DO  100  1=1 .nPS 

IF  (JMEM(l)  .EE.  NCN ) GO  TO  100 
PRINT  90.  I . NCN 

90  FORMAT  ( •0ERRJR....N0DE' . I 5. • 

CALL  EXIT 

100  continue 

yHE.  ORIGIN  OF  THE  NEW  NUMBERING  SYSTEM  IS  LOCATED  AT  EACH  NODE 
IN  TURN  TO  SEARCH  FOR  THE  ONE  WITH  SMALLEST  BANDWIDTH 

DO  160  IK=1,NPS 
DO  120  J=1 . NPS 
JO  I NT ( J ) =0 
120  NEwOT ( J ) =0 
MAx  = 0 
I = I 

NEwOT ( 1 > = 1 X 
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HAS  MORE  THAN'. 15.'  CONN.  NODES') 


■ i'yjZ, 


JO 1 NT ( IK) = 1 
K=1 

1 JO  NEw=NEWOT  < I ) 

K4  = JMEM(NEW ) 

IF  (K4  .EG.  0)  GO  TO  145 
JSUB* (NtwJT ( I ) -1 ) *NCN 
DO  140  JJ=1«K4 
JP=JSUB+JJ 
K5=MEMJT( JP) 

IF  (JCINT(KS)  .GT»  0)  GO  TO  140 
K=K+  1 

NE>|JT(K)=K5 
JOINT (Kb) »K 
ND1FF=1ABS(  I-K.) 

IF  (NDIFF  *GE*  MAXN)  GO  TO  160 
IF  (NDIFF  *GT«  MAX)  MAX=ND I FF 
140  CONTINUE 

IF  (K  .EO.  NPS ) GO  TO  150 
145  1=1+1 

GO  TO  1 30 
150  MAxN*MAx 

DO  155  3=1  « NPS 
155  NJNT( J) =JOlNT < J ) 

160  CONTINUE 
return 
end 
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SUBROUTINE  MPi-0  f I XC  « YC  «NODE  » NN  * NE « NELS  « NPS  > 


f 

I 


C USE  CALCOMP  PLOTTER  TU  PLOT  GENERATED  MESH 

COOED  bV  CApT.  GtWHY  VANKLUREN  OF  AF  F UL » WR  1 GHT-PAT  TEWSON  A(- a 

C 

D I MENS I ON  XC ( NN  ) » YC  ( NN ) * NODE ( NE  » A ) • XP ( 7 ) » YP ( 7 ) 

C DEFINE  PLOT  SIZE  AND  PLACE  PEN  IN  INITIAL  POSITION 

DO  5 L= 1 * NPS 
XC(L)=XC(L)*8. 

5 YC(L)=YC(L)*8. 

CALL  PLOT < 0 • 0 

CALL  PLOT ( 1«0»3*0*“3) 

C DRAW  EACF  ELEMENT  IN  TURN 

DO  10  I=1*NELS 
NPE  = 5 

IF  ( NODE ( 1 « A ) • E U • 0 ) NPE  = A 
NPp=NPE-l 
DO  9 J=  1 « NPP 
XP( J)=XC(N0dEI  I * J > > > 

9 YP( U > = YC  < NODE ( I » J > ) 

XP ( NPE ) = XP  < 1 ) 

YP(NPE)=YP< 1 ) 

XP ( NPE+ 1 ) =YP( UP t+1 1=0.0 
XP<  NPE+2 ) =YP( NPt+2 ) = 1 .0 
10  CALL  LINE <XP»YP» NPE* 1 *0*0) 

WRITE  NODAL  NUMBERS  AND  FINALLY  END  THE  PLOT 
DO  20  L= 1 * NpS 
PL  = FLOAT  <L ) 

20  CALL  NUMBER ( XC ( L> «YC(L) •• 1A*PL«0*-1 ) 

CALL  PLOTE 

return 
END 
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A FORTRAN  listing  of  the  source  deck  for  the  STRANL-II  program  is 
presented  in  the  following  pages.  The  program,  as  presently  configured, 
requires  177gK  words  to  execute. 


PPOGPAM  .'i  A 1 N(  I .\|PUT=ob  .OUTPUT  = 1 3l  .PUNCFi  = bb.  T APLb= INPUT  . 

1 TAPt-6=UU  TPU  T ) 

C 

U jTtrtUY  IpAimScMIC  I-lo*  ANAcYb  I 6 uV  KIMTL  LotMLNT  ml  T HuO 

u-p  I IMU  LLAu  I oUJrt^tJ  * I 1 H T p 1 ANGul Aw  wNj  L.U AOP  I LA  I tPAL  LLc.i’iLi'jTj 

r;  IuPl(l)  = l«  USE.  KtbUL.Tb  OK  PPtVIJUb  CAbh  Aj  SI  APT  I NO  SOLUTION 

C WH 1 t-t  Trio  OTtiLi<  OPTION  lb  IONUkLU 

c i op  i < 2 i = i . ptAi)  i I'j  nun-zlpo  initial  guols 

C 1 uP  T ( 3 ) = 1 . APPL^  LINtAPlZt.0  u>uoNoAPY  CONolTION  ON  liiOPoLINL 

{_  Tr,L  ppcopAivi  P,<c.SLNTLY  UlMLNSlONLO  ALLOWS  T nLSL  MAXIMA 

U 20n  fc.LK  HUNTS'  i Ov)  NOULS.  SO  NouLS  K ok  lACH  TYPL  Oh  UOUNDaPY  LONo. 

C ■ AX  • Full.  uAi'li  a/IJTH  = BA 

C uLVt.LOPLO  ANO  COULO  bY  blT.VLNu  CHAN  OF  LULKHLt  J-HUN  I bV  1 Lll  . AL  A . 

C 

U IMLNj  I 0.4  I I r i-L  l 8 ) • IOP  I <201  «Nuo<  200  «4  > .o  ( b40»64  > . oLP(  S4  0 ) 

O I i4l.No  1 ON  a(  l(jO)  •»(  lbU)  . Pi-iL  l I bO  ) « P/lLP  ( 1 bO  ) « LoF  ( 160  ) «KjinT  ( 160  ) 

UlMLNbloN  N 1 Ob  ( i ) .NIL)  (3»  SO)  . VAF  (.SO)  .API  100) 

LOGICAL  LK(SO) 

LOU  1 VALl-nOL  ( 1. 1 O J ( I ) . I.F  ApF  ) « <Nluo<  2)  .NwAisL  > • < NlOb<  3 ) . NbOoY  ) 

OATA  NLl'l « NPi'1  . Mi'i  .NA/200  . Ib0.64.b0/ 

OATA  Pl/J.  1 4 1 b 9 26/  « GAl'iMA/  I .40/ 

|MKM  = 3 *l  iPM 
I P=2*NA 

CONST  = 0 » S'*  ( GAni-IA-  I .0  ) 
t_Xp  = -GAMMA  / ( GAMMA-  1.0) 

C 

C klAO  TITll.  CoN  r iOl  KLYS.  ANO  PPUdlcM  PAPAmlT LPb 

C 

IoO  PLAO ( b. OOS ) ( r 1 TlL < I ) . I = I * 6 ) 
ik ( t-:uK  < s) ) 2000 »ioi 

10  1 CONTIlTUL 

PC  AO  ( S'  U20  ) ( I Op  r t I ) » 1 = 1 » 20  ) 
plAU  630.  1 To  I V f / 1LST .PnAC » K 1 

PP i NT  Vld.  ( T I T Lb ( I ) . I = 1 ' 6 ) . 2 I Lb T 

Paint  620.  ( I OP  I l I ) . I = I . 20 ) 

I s=0 


pF  T=  I .0 

SoM  AC  =Pi*lAC  * X-  2 

IK  ( I OPT ( 1 ) .LO.  I)  GO  TO  362 
O 

PLAU  ANu  Pp  1 i4  I .-iLbrt  OATA.  bOUNUAPY  i.uULji  pnO  AlkKulL  cLOPL 

C 

KLAU  62  S.  Nb  Lo  . upb  . NbvV  . ( N 1 Ob  ( I ) . 1 — 1 .3) 

(<LAU  62b*  ( ( NoOl  I «J  ) * u=  1 >4  ) » I = 1 .NtLb  ) 

PLAO  640.  ( X ( I ) * V ( 1 ) . I - I .NPb ) 

OO  I 1 0 1 = 1.3 
Nb=NlDb( I ) 

110  PLAO  62b.  <NlD( I »J> » J=I *NS) 

PLAO  640.  ( VAF ( I ) . I = I .NbOOY ) 


80 


1^0  kC  Au  82b*  (NJnT  t 1 ) * 1 = 1 .NPb) 


C 

C 

C 

C 


C 

C 

C 


It-  1 UPT  ( 3 ) = 1 » APPLY  LINEAklZLU  UUUNDAkY  LUNDlTlUN  UN  lHoNDLINL 
OlHLkWlSt-  APPLY  NUNL INEAk  bUUNO AkY  LUNJ  I T 1 ON  UN  A I Hi  U i L SC'HIAli. 


11  ( 1 UPT  < 3 ) .IML.  1 > W TU  lib 

DU  lib  ^1  .NODDY 
1 =NlD(3. J > 
lib  Y ( I ) =0*0 
1 lo  LUNl INUE 

200  Pk  1 NT  920*  NtLb.NPb.NBW 


PklNT  930 
DO  220  N = I ,NLLS 

22 0 Pk 1 NT  825*  P » ( NUD (N«J)iJ=l*4) 

Pk 1 NT  93b 
UO  230  I = I .NPb 

230  Pk I NT  940 . I . NJNT < l ) . X ( I ) . Y ( I > 
Pk  I NT  9b  1 < ( |\|1D(  1 « 1 ) . 1=  1 »NP  AMP  ) 

Pk  1 NT  9b2  . (NlD<  2*  l ) . 1=  1 .NWAKE) 
Pk  I NT  9b  J . (NIDI  3. 1) . 1 = I .NbODY) 
Pk  1 NT  9bb . { VAF ( I ) . 1= 1 .NBODY > 


HtDEF INl  PILSH  DATA*  ETC*  U2  1 NO  NEW  NODAL  NUMuEklNG  SYSTEM 

DU  238  N= 1 * Nt  LS 
Du  238  1=1.4 

IP  (NODIN. I ) *E0.  0)  GO  TO  23H 

KK=.NOD<N»  1 ) 

NUD<N. I )=NJNT ( KM 

238  CUNTINUL 

DO  239  1=1.3 
1 S=NIOS< I ) 

DO  239  J=l  . lb 
KK=NID< I « J) 

239  N ID  < I * J)  =NJnT  ( Kk  ) 

DO  243  I = I . NFS 
kMLPl  I ) =X<  I ) 

24  3 k.vin  I >=Y(  I ) 

DO  244  11  = 1 ,NPS 
I =NJNT< 1 l ) 

X ( 1 ) =HMLP(  I I ) 

244  Y ( I ) = HML (11) 


C 

C HALF  BANDWIDTH  AND  NUMBEW  OF  EOUATIUNb 

C 

nhbw=nbw/2 

NEQ=3*NPS 


C kt.AU  NONZEkO  INITIAL  GUESS  Ok  PkOCEED  WITH  ZEkO  SOLUTION 

C 

IF  (I  OPT (2)  • NE • I)  GO  TO  250 
kEAD  840.  ( S < I »N8W ) . I = I * NEU ) 

GU  TO  29b 

2b0  DU  260  1= l ♦ NtU 
260  SLp ( I ) =0.0 


C 

C ITEkATIONS  STaKT  Ht-kE  AND  CHECKED  IF  ITGlV  IS  EXCEEDED.  IF  SO. 

C Pk I NT  FAIL  TO  ClUVeHGE  AND  PkOCEED  TO  NEXT  CASE.  OTHLkWlSE. 
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- - - •— - 


■ - _ . mm.  m . 


CUN  I INUL  ro  | I LKATL 


c 

c 

c 


20b  Ik'(  ;s  = lRLb  + 1 

It  ( IKL3  • Oj  • lTulV)  GO  TO  600 
I-UHMULAIl  SYS  I Li"l  of-  ALGEBRA  I C EUUATlONS 
UO  266  1 = 1*  Nto 

Do  266  u=l  ,nBW 
260  S ( I . J ) =u. 0 

call  NEWN(  SGli«iAC  » FIRM  *NCM  • NEU*  NbW  « NEM  « NELS  * NOD  * aLP«  S»  COE  ,,\|PM,x  . Y ) 

IMPOSE  u.  C.  hOK  KARFlfcLu.  LINE  OF  SYMMETRY,  AND  ON  AIRFOIL 

OO  274  1=1 ,NFARF 

1E=J*N1U< 1,| ) -3 

OO  272  N = l,3 

It-=1E+1 

DO  270  K=  1 , NBtv 
£<0  S(  IE  ,K  ) =0.  0 
B7 c'  S<  I E » NHLSW  ) = | ,0 
2 74  CONTINUE 

OO  280  I =1  * Nw Ak  c 
I L = 3*N 1U  < 2 , I ) 

OO  278  k=  1 ,nB/J 
<^?8  S < IE  *K  ) =0»  0 
280  S(  IE ,NHbW  ) = l .0 

OO  280  U= 1 , nBOOY 
1 E = 3*N 10(3 , j ) 

DO  282  (2=  1 , nUW 

2el2  S(  IE  ,K  ) =0. 0 

S(  IE,NHUw)  = l .0 

IF  ( I OPT  ( 3 ) «i\i£,  l)  S < l E *NHbw- 1 ) =-VAF ( J > 

26b  b ( | E * Nb*l ) = VAF  ( J ) 

CALL  bANDEU  SOLVER  TO  SOLVE  THE  SYSTEM  OF  LOCATIONS 
CALL  UNOEO  < s » mUi'i  *NCM  , NEU ,NHBw  ) 

PRINT  COMPUTED  KtSuLTS 

2Vb  PRINT  970,  RK  AC  » RFT 
PRINT  97b,  IKES 

00  Job  I = 1 , nPS 
J=NJNT ( I ) 

1 I = J*NJNT ( I ) 

POT=S<  I 1-2. NBW) 

UPT=S<  I 1-1  ,NB*/  ) 

V = S(  I I , NbW ) 

U=uPT+ 1 .0 
USQ=U*U+V*V 

AbQ  = CONST'*-  ( 1 • 0-USU  ) + 1 ,0/SQMAC 
RML( J)=SORT( OS0/ASO) 

PRAT  I 0=1  1.0+CONST  *RML <J)**2)**EXP 
CP=— 2.*UPT 

COF ( J ) = SUMAC  * ( 1 • 0+2,4+UPT ) 

DELM=RML ( J ) — RmLP ( J ) 

80b  PRINT  978,  I ,P0T,U,V,COF(J),RML(O),PRATlO,CP,OELM 

oi splay  pressure  coefficient  cp 

PRINT  96b 
I STOP=0 
I PNT  = 0 
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fW 


DO  320  J=  1 .NbODY 

IF  (J  .t<J.  NbODY)  1ST  OP  = 1 

1 =NID( 3. J) 

UPT=S(3*l-l,NbW) 

V=S ( 3* 1 .NB* ) 

CP=-2.*UPT 

820  CALL  FPLUT ( I A. IPnT.AR.LR. ISTOP.  1 . 1 .X(  I ) .CP) 

IF  < lktS  .LT.  2)  CO  TO  382 
c 

C CHECK  CONVERULlNEt.  IF  SO.  PUNCH  CONVERGEu  SOLUTION  ANU 

C PKOCttU  Tu  NLAT  CASE.  OTHERWISE*  UPDATE  SOLUTION  AND  CONTINUE 

C TU  1TEWATE 

C 

DO  340  1 = 1 .NPS 

PCTE=1  .0-Hi>ilP(  I )/HML(  I > 

IF  (ABS(PCTF)  .LT.  ZTEST ) GO  TO  340 
HFT=F1 
CO  TO  JB2 
340  CONTINUt 

PUNCH  040.  (SU.NBW).  1 = 1. NED) 

GO  TO  100 

303  UO  30b  1 = 1 .NPS 
38b  PMLPI I )=WML< I ) 

RFTC= 1 • O-RFT 
DO  300  1 = 1 . NEo 

390  SLp  ( 1 ) =RF  T*S  ( I » i 'ItJW  ) +RF  TC*SLP  ( I ) 

GO  TO  265 
600  PU I NT  980 
GO  10  10U 

805  FORMAT  < 7A 1 0 , A 2 ) 

Ui.0  FORMAT  (4012) 

8c 5 FORMAT  (1615) 

BsO  FORMAT  (I5.3F10.0) 

840  FORMAT  (IPbElO.ii 

910  FORMAT ( 1H1 «2A » 7A10* A2// • CONVERGENCE  LImIT  =».F6.4//) 

92C  FORMAT  (1H0.-IM0.  OF  ELEMENTS=*  * 14.  • NO.  OF  N0DES-'.I4, 

1 • FULL  DANDWIDTH= • , 14//) 

930  FORMAT  ( ' OELt • NO.  AND  ELEMENT  NODES'/) 

935  FORMAT  ('OULO  NODE'.'  NEW  NODE*.  6X . • X ( I ) • . 1 2X . • Y ( I ) • / ) 

940  FORMAT  ( 16. I 10  » 2E 1 5.5 ) 

951  FORMAT  ('ONODLS  AT  FARF I ELD • / < 20  1 5 ) ) 

982  FORMAT  ('ONODLS  ON  LINE  OF  SYMME TRY • / ( 20 1 b ) ) 

953  FORMAT  ('ONODCS  ON  THE  A I RFO I L • / ( 20  I 3 ) ) 

98u  Format  ('osloPE  Along  nodes  on  a irfo i l • / ( 8l i 5.5 ) ) 

9/0  Format  <1H1,  -MACH  NUmBER=  • . F 6 . 3 « bX  . • Rt  LAX.  FaOToR  = '.Fc>.4/> 
9/5  FORMAT  (1H0.4A.-NO.  OF  I TERAT I ON= • * I 4/ 

1 1 HO  , 7 X .-NODE  • iBX.  • PFI  IT'  . 1 IX.  • UCOM  • . 1 1 X . • VCUM  • , 

2 I2X.  'COF - , 1 IX.  • LMAC • . 1 IX.  • P/PU • . 1 3X . • CP  * . 1 1 X . • DcLM • / ) 
978  t-ORMAT  < 1 10,lPbci5.4> 

98C  FORMAT  ('OFAlL  TO  CONVERGE  IN  SPLCIFILD  NO.  OF  ITERATIONS') 
985  FORMAT  ( 1 FI  1 ) 

2000  STOP 
END 
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bUBWOUT I NE  NE WK 1 SUMAC  « NHM « NCM * N£U * NBW  * NEM . NELS  * NOD  * 3LP  * S * 

1 COEF,NPM.X»Y  ) 

OLnERATL  SYSTEM  MATRIX  BY  AS5EMBLL1NG  CONTRIBUTIONS  f-NOM 
ALL  THE  ELEMENTS 

SUBROUTINE  EMuT  CALLED  TO  GENERATE  ELEMENl  MATRIX 

DIMENSION  COEK  1 NPM) «X<NPM> «Y<NPM> . Xu< A ) * YU<  4 ) • PMC  12)*  BB  <12*12) 
D 1MENS10N  NCU  l NEM*  4 ) «5<NRM*NCM)  • i»LP<NRM  ) 

NHBW=NBW/2 
DO  480  Ns  1 « NELS 
11  = 1 

IF  < NOD ( N * 4 ) ) 402*402*404 

402  NPEL*3 
NTRS* 1 
GO  TO  410 
404  NPEL*4 
NTRS *4 
12  = 0 

DO  408  1-1*4 
Nl=NOD<N«  1 ) 

408  IF  CCOEFCNl)  «GT«  1*00)  12=12+1 

IF  <12  .EG.  0 ) NTRS=2 
IF  < 12  .EG.  4)  11=3 

410  DO  425  1=1 .NPEL 

Nl=NOD<N, 1 ) 

XU< 1 ) =X  < N 1 ) 

YU< 1 )=Y(N1 ) 

DO  425  J= 1 *3 
1 S=3*  < N 1 — 1 )+J 
1E=3*< 1-1 )+J 
425  PM< IE ) = SLP< 1 S ) 

CALL  EMUT  < XQ  * YG  i PM*  SQMAC *88* NTRS ) 

DO  450  1=11 .NPEL 

NR=3*<NOD<  N , I ) -1 ) 

I E = 3*  < 1-1  ) 

DO  450  11=1*3 

NR=NR+1 

1E= 1E+1 

DO  450  G= 1 «NPeL 
NC=3*(NOD<N, J ) - 1 ) -NR+NHBW 
JE  = 3*<  J-l  ) 

DO  450  JJ  * 1 , 3 

NC=NC+1 

JE= JE+1 

450  S<NR«NC) = S< nR»NC )+BB< IE *0E ) 

4B0  CONTINUE 
RETURN 
END 
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.1,.^ 


c 

c 

c 

c 


SUBROUTINE  eMUT  < XU. YU *PMU » SUMAC *EU « NTRS ) 


GENERATt.  MATRIX  FUR  A UUADR1 LATERAL  OR  TRIANGLE 

SUBROUTINE  EM  f C CALLEU  TO  GENERATE  MATRIX  FOR  A bASlC  TR 1 ANGLl 


iOO 


102 


105 


DIMENSION  EQ< 12  » 12) »ET(9»9>  «XU(4 ) * YU(4) *XT (3) *YT( 3) »MP( 3*4 ) 
DIMENSION  PMU( 1 2) *PMT (9) 

DATA  MP/1 .2, 3 *3 .4* 1 .2.3*4. 4. 1 *2/ 

FTOR  = 1 .0 

IF  (NTRS  .EQ.  4)  FT0R*0*5 

DO  100  1=1.12 
DO  100  J=1 «i2 
EU( 1 , J)=0.0 
DO  150  11=1, NTRS 
DO  105  1=1,3 
N 1 = MP ( I , I I ) 

IT=3«< 1-1 ) 

I 0 = 3* ( N 1 — 1 ) 

DO  102  J=  1 * 3 
IT=IT+1 
lu= 10+1 

PMT (IT) =PMU (10) 

XT ( I )=XU(NI  ) 

YT ( 1 ) = YU ( N I ) 

CALL  EMTC(ET»XT,YT.PMT.SQMAC> 

DO  130  K= 1 ,3 
NR  = 3« (MP(K , I I) -1 ) 

1E=3# (K-l ) 

DO  130  KK* I , 3 

NR=NR+1 

IE= IE+1 

DO  130  L=I ,3 

NL  = 3* ( MP( L , 1 11  — 11 

JE=3*(L-l ) 

DO  130  LL* 1,3 

NC=NC+1 

JE=JE+1 

130  EU(NR,NL)sEQ(NK,NC)+ET( 1 E » JE ) *F  TOR 
150  CONTINUE 
RETURN 
END 
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subroutine  emtci A.XL. YL.PEL.SUMAC) 


c 

C EVALUATt  ELEMtN  T MATRIX  FOR  A TRIANGLE  BY  GAUSSIAN  GUADRATURE 

C SUBROUTINE  dE«V  CALLED  TO  EVALUATE  SHAPE  FUNCTION  DERIVATIVES 

C AT  THE  GAUSSIAN  POINTS 

C 

DIMENSION  A ( 9 * 9 ) *P I 9)  .0t9>«NP<5>»B<3)*Ct3>**L‘3>«YL(3>*Sl3>* 

1 DnX ( 9 ) »DNXX ( 9 ) «DNYY l 9 > .PEL (9) 

DIMENSION  EINII3.7).  WT ( 7 ) 

DATA  LMAX/ 7/ > H T / 0 • 225 . 3*0 • 1 32394 1 b . 3*0 • 1 259J9 1 ti/ 

DATA  EINT/J#0»  33333333.  0 • 05961387.3*0* A 70 1 4206 » 0. 05961 5B7 . 

1 3*0.47014206  » 0 • 0596  1 56 7 » 0 • 79743699  « 3*0  • 10126651  . 

2 0. 79742699*3*0. 1 0 1 2665 1 . 0 • 79742699/ 

DATA  NP/ l*2«3»l*2/» GAMMA/ 1.40/ 

DO  2 1=1.9 
DO  2 J= 1.9 
2 A( I . J)^0. 

DO  4 1=1.3 
J=NP< 1+1 ) 

K=NP< 1+2) 

IU  I )*YL(  J>“YL(ia 
4 C(  I ) = XL  IK) —XL  I J ) 

ARE A = 0»  5*  I B ( 2 ) *C ( 3 ) — B ( 3 ) *C  < 2 ) ) 

CST 1*1 .O-SGmAC 
C5T2=SOMAC* ( 1 . O+GAMMA ) 

DO  100  L= 1 .LMAX 
DO  10  1=1.3 
10  S(  I ) *E 1 NT  < I ,L) 

call  derv < area. b.c.s. dnx. dnxx . dnyy > 

u = 0. 

ux=o. 

DO  30  1=1.9 
U=U+DNX( I )*PEL( 1 ) 

30  UX=UX+DNXX< 1 )*PEL< 1 ) 

ALPHA=CST1-cST2*U 
DO  40  I = l « 9 

P I I ) =ALPHA*dNXX ( I ) +DN YY  < 1 ) 

40  0(  I ) =P< I >-CsT2*UX*DNX ( l > 

XEIGHT=WT(L)*AREA 
uO  60  1=1.9 
CST=WEIGHT*Q( I ) 

DO  60  J= 1 . 9 

60  A( I , J)=A( I , J)+CST*P< J> 

1O0  CONl INUe 
RETURN 
END 
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SUBWOUTlN  DtKV  ( AREA . b . C « S «DNX « DNXX »DNYY  ) 


t YAL.UA  1 E 3HAPL  f-UNdlON  DEMIVAllVEB  AT  OAUBBlAN  POINT 

OlMENblON  t J < 3 ) » C ( 3 ) « S < 3 ) «DNX<9)  «DNXX(9)  «UNYY<9>  *NR<5> 
OAT  A NR/I *2, 3 *1 *2/ 

TW0A=2.*AREA 
TWOASU=TWOA**2 
DO  200  I = I • 3 
J = NP<  I + I ) 

K.  = NP<  1+2) 

Sl=S< 1 ) 

SJ=S<  J) 

SK=S(K) 
fc?  1 =B  ( 1 ) 
bJsB( J) 

BK=b(K) 

C1=C(  1 ) 

CJ=C<  J) 

CK=C(K) 

i>  1 SU  = b I *S  1 

BlSCJ=Bl*Bi 

ElSQ=Cl*Cl 

ALFA  = 0»b*<  CK-CJ ) 

BETA»0.b*  < BJ-BK.  ) 

HX=BI*SJ*5K+B J*SK*SI+BK*bI*SJ 
HXx=2«*<SI *BJ*8N+bJ*BK*BI+SK*BI*bJ) 

HYy=2.*<S1 *CJ*CN+SU*CK*C1+SK*CI*CJ) 

CSS=6.*S1 * < 1 .-S I ) 

CS=6.*< 1 .-2.*SI ) 

-=3*1-2 

DNx i L ) =B I *CSS+  2 . *HX 
ONxX ( L ) =B I SQ*CS  + 2 • *HXX 
DNyY<L)=CI  SQ*CS  + 2»*HYY 
CS  = CK*SU-CJ*5K. 

L=E  + 1 

DNx  < L ) =2.*Bl*Sl*tS  + TWOA*bISO+ALFA*HX 
DNxX ( L > =2. *B I SQ*  CS  + 4 . *B I *TWOA*S I +ALFA*HXX 
DNYY  < L ) =2. *C I SQ*CS+ALFA*HYY 
Bb=BJ*SK- 3K*S J 
L-=L+  1 

DNx ( L ) =2 • *B I *S I *BS+BETA*HX 
DNXX < L ) =2. *B I SU*BS+BETA*HXX 
200  DNyYIL) =2.*C ISU*BS+4.*Cl*TWOA*SI+bETA*HYY 
DO  J00  1 = i 19 
ONX ( I >=UNX< I ) / TWUA 
DNXX ( I )=DNXX( I ) /TtfOASQ 
300  DNYY ( 1 ) =ONYY ( I > / TWOASQ 
RETURN 
END 


SUBROUTINE  BNOEUI A.NRMAX.NCMAX.N.  1 TERM ) 


C 

C EQUATION  SOLVER  FOR  BANDED  NON-S YMMETR I C SYSTEM  OF  EQUATIONS 

C SOLUTION  STORED  IN  THE  LAST  COLUMN  A(  1.2*1  TERM) 

C 

dimension  A(NKMax.ncmax) 

CERO= 1 »E-6 
PARE  = CERO*#2 
Nt)ND  = 2*  1 TERM 
NBM  = NbND  - 1 

C BEGINS  ELIMINATION  OF  THE  LOWER  LEFT 

DO  1 OOO  I ~ 1 , N 

IF  ( ABS<  U I . I TERM) ) .LT.  CERO ) GO  TO  410 
GO  TO  430 

410  IF  < ABS<A( 1 , IJERM) ) .LT.  PARE)  GO  TO  1600 
PRINT  420.  A ( I • I TERM)  . 1 

420  FORMAT  (•  WARDING*  1 LL-CONO I T I OREO  A-MATRlX.  A**,EI6.6»*  I=*.l4) 

4j0  JLAST  = M1N0I  1+ I TERM- 1 . N) 

L = I TERM  + 1 
DO  500  J=1 » JLAST 
E = L - 1 

IF  < ABS ( A < j . L ) ) .LT.  PARE)  GO  TO  500 
B r A ( J ,L  ) 

DO  450  <=L.  NBND 
450  A(J.K)  = AIj.K)  / B 

IF  <1  .EQ.  N)  GO  TO  1200 
500  CONTINUE 
0 = 0 

OFIRST  =1+1 

IF  < JLAST  ,Lt.  I)  GO  TO  1000 
DO  900  J=  JFIKST.  JLAST 
L = L + 1 

IF  < ABS< A( j, ITERM-L) ) .LT.  PARE)  GO  TO  900 
DU  600  K=ITeRM,  NBM 
t>O0  A(J.K-L)  = A ( J — L *K)  - A(J.K-L) 

A (J. NBND)  = A( J-L.NBND)  - A(J.NBND) 

IF  (I  .GE.  N-ITERM+l)  GO  TO  900 
DO  BOO  K=  1 « L 

800  A(J.NBND-K)  * -A(J.NBNO-K) 

900  CONTINUE 
1000  CONTINUE 
C oACK-SUbST I T U T I UN 

1200  L = I TERM  - 1 
DO  1500  1=2,  N 
DO  1500  J= 1 , L 

IF  (N+l-l+J  .GT.  N)  GO  TO  1500 

A ( N+ 1 — I .NBND ) = AIN+l -l .NBND ) - A < N+ 1 - 1 + J , NBND ) *A < N+ 1 - I . I TERM  + J ) 
1500  CONTINUE 
RE  TURN 

C PRINT  THE  ENTIRE  matrix  IF  ZERO  ON  main  DIAGONAL 

1600  PRINT  1601, 

1601  FORMAT  !'  COMPUTATION  STOPPED  IN  BNDEG  BECAUSE  ZERO  APPEARED  ON 

1 MA I N DIAGONAL.  THE  MATRIX  FOLLOWS.*) 

DO  1602  1=1,  N 

1602  PRINT  1603.  (A(I.J),  J=l,  NBND) 

1603  FORMAT  I 10E12.4) 

STOP 

END 
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SUBROUTINE  FPlOT  (Ml  . 1 PNT  « AR  * l_R  * 1 STOP .NC * NCMAX . V 1 . V2 ) 

POINTER  PLOT  FOK  DISPLAYING  COMPUTtD  RESULTS 

LOGICAL  LmI  2 ) .t-N(  120)  *LP<4  ) »LX(  4 ) *LR(  1 > .LC(  20) 

INTEGER  R1  (2)  .ST12) .SG<2) 

dimension  AR(  mi  ) «N(  120)  • WHO  t 30 ) «/< IS) • 1 OFE  l 2 ) « l GP  ( 2 > * sF  1 l 2 ) « Sf- 2 < 2 ) 
1 » I D 1 2 > » I 5u ( 2 ) • I l ( 4 ) . 1 P ( 4 ) » 1 M ( 4 ) • 1 X ( 4 ) « IR(b02 ) » 1 A( 23U) 

Euu  1 V ALLNCL  ( IP  ( 1 ) *LP  ( 1 > ) • ( 1M(  1 > »LM  ( 1 ) ) • < LN  ( 1 ) <N  ( 1 ) *RHO<  1 > ) » 

1 ( 1X<  1 ) *LX( 1 ) ) . ( N120.R1 < 1 ) > • (N5B.RI  ( 2) ) 

DATA  (/(  1 ) • I = 1 • it)  >/  1 ■ • 1 • 25  . 1 * 5 * 1 ■ 7b  « 2.  « 2 • 5 • 3.  .3.5*4. »4.5*5.*b.« 

I 7»  »8*  «y •/ 

DATA  < 1 SP ( 1 ) . 1 = 1 « 2 ) » ( I X ( 1 ) ♦ 1 = 1 • 4 ) • ( I P ( 1 ) * 1 = 1 . 4 ) . ( 1 M ( 1)  • 1 = 1*4)/ 

1 10*5*'  »lHX«4*lH+.4*0/ 

DATA  1D( 1 ) . 1 D 1 2 ) ♦ 1SG(  1 ) . I SGI  2)/ 120*58* 1.-1/ 

DATA  LNC 1 > . 1 1 »LC( 1 )/l .80*50.0. 1 . 1H0/ 

1F( 1 STOP ) 173.172.172 

172  J = I PNT +2 

1K( J.GT.M1 )G0  TU  173 
ipnT  = j 
AR( J-l )=V1 
AR( J)=V2 
DO  10  K=1 .4 
lO  1 M( K ) = NC 
J=J/2 

LR( J) *LM(2 ) 

173  1F(NC. LT. NCMAX  1KETURN 

1F(  ISTOP.EQ.O IRE  TURN 
DO  171  1=1.? 

R 1 ( I ) « 1D( 1 ) 

SG(  1 ) = 1 SG( 1 ) 

IF(  1 1 ( 1 ) .NE.O 1R 1 l 1 > = 1 1 ( 1 ) 

IF  ( 1 I ( 1+2) . EQ. 1 ) SG<  I >=-SG<  1 ) 

171  ST(  1 ) = (W1 < I ) + l-SG(  1 >*CR1 ( 1 )-l ))/2 
1F(R1(  1 ).GT.120)«1  ( 1 > = 120 
1F(R1(2).GT.238)R1 (2) =238 
N2o  = N 1 e 3/6 
N 1 20  = 6*  N20 
DO  3 1=1.2 
AMI  Ns AR ( 1 ) 

AMAX  = AR ( I ) 

DC  7 J=l.  1 PnT  . 2 
A A 1 = AR  ( J ) 

1F< AMIN.GT.AA1 )AM1N=AA1 
7 1F< AMAX.LT. aAI )AMAX=AA1 

IF  (ASS(AMIn)  .LT.  .000001)  AMIN  = 0.0 
IF  (ABS(AMAx)  .LT*  .000001)  AMAX  = 0*0 
R=AMAX-AM1N 

IF  ( ABS(R  ) »LT.  1 • t-9. AND. AbS( AMAX) .LT. 1 .E-9)  R=1»0 
IF  ( ABS(R)  »LT.  1 . t-9. AND. AMAX.LT.O. ) R = -AMAX 

IF  ( AbS(R)  »LT  • 1 . t-9. AND. AMAX.GT.O. ) R=  AMAX 

DO  22  J=  1 . 15 

B = ALOGl 0 ( R / ( R I ( I ) -2 ) /Z ( J ) ) 

M = b 

IF ( B.LT *0 ) M=M— 1 
C=Z( J)*10.**( M+l  ) 

B = AMIN/2(  J)/10.**(M+1  > 

1 1 =d 
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IF(B.LT.O)  I 1*11-1 

IF ( ( R I ( I 1-2+11 )*C-AM1N)  18*19.19 

18  C=10.*C 

19  1F( J.EQ.l >SMlN=C 

22  1F(C.LT.SM1N)SMIN=C 

SF  i ( | ) = ( 1./SMIN)*SG< 1 1 
B=AMIN/SM1N 
M = 0 

1F(B.LT»0)MsM-1 
bF  2 ( I )=bT<  I )-M*SG(  I 1 
WHO l I )=SF2< l ) +0 • 5 
M=sF2< 1 ) 

DO  25  J=1 » 10 

IF  ( M— J—  ( <M-J)/  Ii>P(  I ) )*1SP<  1 > 125.3*25 
25  CONTINUE 
3 I OFF  < I 1 = J 

DO  101  1*1. N58 
10  1 IA(  I 1 =0 

DO  102  J= 1 , IPNT  *2 
I W ( J 1 =SF 1 ( 1 i *AR 1 J 1 +RHO ( 1 1 
1 T = SF 1 <21*AR( J+l ) +RHO ( 2 ) 

1 F ( J • NE • 1 1 GO  TO  109 
IK( IPNT+2 1 =2 
I«( J+l 1 *0 

108  13= IT 

GO  TO  102 

109  1F(  IT-131  104*105* 105 

104  1R< J+I 1=IR< lPNT+21 
1 R ( I PNT +2 1 * J+ 1 

GO  TO  108 

105  1 * I T+ I 

106  1=1-1 

1 1= 1A( 1 1 

IF( 1 1 1 106. 106.107 
107  1W( J+l 1=IR< 1 1 ) 

1W( 1 1 1=J+1 
102  I A ( I T 1 = J+ 1 
LAST= IPNT+2 
J J= 1 OFF ( 2 1 
LZH=SF2( 1 1 
LZv=SF2<2! 

DO  100  1*1. N58 
DO  40  J=l.Nl20 
40  N( Jl = 1H 

IF ( ( 1-1 1*( 1-N58 1)140.151.140 
151  DO  141  J= 1 «Nl 20 
14  1 N < J ) = 1X(  1 ) 

140  LN< 1 ) =LX< 1 ) 

LN( N 1 20 ) = -X ( 1 ) 

1F(LZH.LE*O.OK.LZH.GT.NI20)GO  TO  131 
LN( LZH) =LP<  4 ) 

131  NB=1 

1 F ( 1 .NE.JJ  >G0  TO  35 
Nb=2 
J J= JJ+5 
1 3*1  OFF ( 1 ) 

DO  32  J* 1 3 « Nl 20  *10 
32  LN(J)=LP(4> 
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IF ( I .NE.LZV )G0  TO  35 

00  135  3=1 .Nl20 
135  N( J>  = 1P(  1 > 

35  13= 1 A< 1 ) 

I F ( I 3.EO.0 )G0  TO  121 

120  LAST*IW(LAST) 

12= IR1LAST-1  ) 

I l=LAST/2 

LM<  2) =LR<  I 1 ) 

1 MM  = I M ( 2 ) 

LN( I2)=LC< IMM) 

1F(LAST.NE.|3)G0  TO  120 

121  GO  TO  (38.41 ) »NB 

38  WRlTE<6*39)  (N<J)  #J=1»N120) 

39  FORMAT  < 1 I H .120AJ1 

GO  TO  100 

41  AA 1=1 

VALUE=<AA1-SF2<2 ) l/SFl ( 2 > 

U/Rl  TE<  6.42  1 VALUE  * ( N<  3 ) « J = 1 .N I 20 ) 

42  FORMATt  1H  , 1 P 1 1 0 • 3 . 1 20 A 1 ) 

100  CONTlNUt 

1 3= 1 OFF ( I ) 

3 = 0 

DO  49  1=13. Nl 20 *10 
3 = J+  1 
A A=  I 

49  RHO(3)=(AA-SF2( 1 ) J/SFl ( 1 ) 

I F ( 1 OFF ( 1 1-5)62*62.63 

62  */R  1 TE  < 6*50  ) ( RHO  ( I > « I = 1 « 3 ' 

50  FORMAT< 10X, 12( 1PE10.3) ) 

RETURN 

63  IF(J.GE. 12)0=11 

«/Rl  TE  (6*64  1 ( RHO  ( I ) « 1 = 1 *31 

64  FORMATt 16X, l 1 ( 1PE10. 3) > 

RETURN 

END 
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